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Abstract. This paper presents a complete theoretical framework for plasma 
turbulence and transport in tokamak plasmas. The fundamental scale separations 
present in plasma turbulence are codified as an asymptotic expansion in the ratio 
e = pi I a of the gyroradius to the equilibrium scale length. Proceeding order-by- 
ordcr in this expansion, a framework for plasma tiirbulcncc is developed. It comprises 
an instantaneous equilibrium, the fluctuations driven by gradients in the equilibrium 
quantities, and the transport-timescale evolution of mean profiles of these quantities 
driven by the interplay between the equilibrium and the fluctuations. The equilibrium 
distribution functions are local Maxwellians with each flux surface rotating toroidally 
as a rigid body. The magnetic equillibrium is obtained from the generalized Grad- 
Shafranov equation for a rotating plasma, determining the magnetic flux function from 
the mean pressure and velocity profiles of the plasma. The slow (resistive-timescale) 
evolution of the magnetic field is given by an evolution equation for the safety factor 
q. Large-scale deviations of the distribution function from a Maxwellian are given 
by neoclassical theory. The fluctuations are determined by the high-flow gyrokinetic 
equation, from which we derive the governing principle for gyrokinetic turbulence in 
tokamaks: the conservation and local cascade of the free energy of the fluctuations 
(i.e., there is no turbulence spreading). Transport equations for the evolution of the 
mean density, temperatiire and flow velocity profiles are derived. These transport 
equations show how the neoclassical and fluctuating corrections to the equilibrium 
Maxwellian act back upon the mean profiles through fluxes and heating. The energy 
and entropy conservation laws for the mean profiles are derived from the transport 
equations. Total energy, thermal, kinetic, and magnetic, is conserved and there is 
no net turbulent heating. Entropy is produced by the action of fluxes flattening 
gradients, Ohmic heating, and the equilibration of interspecies temperature differences. 
This equilibration is found to include both turbulent and coUisional contributions. 
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Finally, this framework is condensed, in the low-Mach- number limit, to a concise set 
of equations suitable for numerical implementation. 



Multiscale Gyrokinetics for Rotating Tokamak Plasmas 



3 



1. Introduction 

Plasma turbulence in fusion devices is a fundamentally multiscale problem both in 
space and in time. The turbulent fluctuations driven by background gradients typically 
occur at scales associated with the ion Larmor radius, pi (or smaller), whilst the mean 
temperature, density, and bulk velocity profiles vary smoothly over the system scale 
(in tokamaks, the minor radius a). Similarly, the fluctuation frequency u is much 
larger than the rate at which the mean profiles evolve, ~ t^^, where te is the energy 
confinement time. Table 1 gives approximate values for these space and time scales in 
some large tokamaks - it is manifest that the separation of scales is very strong in such 
plasmas. Because of this scale separation, it is possible to average over the time and 
space scales associated with the turbulent fluctuations and consider mean fields that 
slowly evolve due to turbulent and coUisional transport, with the turbulence in turn 
driven by gradients in the same mean fields [1, 2, 3, 4, 5, 6]. 

There are also scale separations associated with the turbulence itself. Turbulence 
in a strongly magnetized plasma occurs at frequencies u which are much smaller than 
the cyclotron frequency of the ions, Qi (see Table 1). The turbulence is also strongly 
anisotropic, viz., correlation lengths along the mean magnetic field are much longer than 
correlation lengths across the field; particles can stream rapidly along field lines but only 
drift slowly across them. These two properties of the turbulence are the foundation of 
the gyrokinetic theory [7, 8, 9, 10], in which the fast cyclotron time scales are averaged 
out and the full 6D kinetics reduced to a simpler 5D formulation, the kinetics of charged 
rings. 

In this paper, we unify this hierarchy of timescales and spatial scales in one 
formulation. We use the physical scale separations inherent in plasma dynamics to 
determine how the mean fields influence the evolution of the small-scale turbulence, 
and how the turbulent fluctuations react back upon the mean fields. To quantify this 
back reaction, we derive the transport equations, in which it becomes manifest by what 
physical mechanisms the turbulence can affect the mean distribution functions and fields. 

The structure of this paper is as follows. We start by recapitulating the fundamental 
equations of plasma physics in Section 2 and splitting all quantities into mean and 
fluctuating parts in Section 3. We then impose order on our multiple small parameters 
and scale separations by introducing the fundamental gyrokinetic ordering in Section 3.4. 
This allows us to formulate the entire problem as a systematic expansion in the small 
ratio e = pi/a. In subsequent sections, we proceed to expand the equations of Section 2 
order by order in e, zeroth order in Section 4, first in Section 6, second in Section 7 
and finally, third (transport order) in Section 8; pausing in Section 5 to introduce 
rotating gyrokinetic variables, which are convenient for handling the turbulent kinetics 
in toroidally rotating plasmas. 

Sections 2 to 8 are, necessarily, quite technical. Readers who believe themselves 
to be already familiar with the formalism presented in these sections may wish to skip 
directly to Section 9. Sections 9 and 10 attempt to explain the physical content of 
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multiscale gyrokinetics by combining all our earlier results and examining how energy 
and entropy flow between the various constituent parts of our multiscale system. In 
Section 9.1, we prove that our transport equations conserve the total energy, and that 
the fluctuations can do no net work on the mean flelds. Explaining this result leads 
us in Section 9.2 to consider the balance of free energy of the fluctuations, as this is 
the fundamental conserved quantity of kinetic turbulence [11, 12, 13, 14]. The equation 
derived for the conservation of free energy clearly demonstrates a local (to a given flux 
surface) turbulent cascade that takes the energy injected by instabilities and dissipates 
it via collisions. In Section 10, we link the free-energy cascade and the mean-fleld 
transport through the evolution of the mean entropy of the system. 



Parameter 


JET 


D-IIID 


K-STAR 


TFTR 


ITER 

(projected) 


B, T 


3.5 


2 


3.5 


5-6 


5 


rig, cm^^ 


5 X 10^3 


5 X 10^3 


5 X 10^3 




10^4 


Ti, keV 


5-15 


7 


5-10 


5-32 


25 


u, km/s 


500 


200 


0.1-100 


100 


50 


M = u/v^Yy, 


0.3-0.5 


0.35 


0.01-0.2 


0.1-0.2 


0.05 


Aoe, cm 


1.3 X 10-2 


8.8 X 10-3 


8.2 X 10-3 


4 X 10-3 


1.2 X 10-2 


Pi, cm 


0.051 


0.06 


0.045 


0.032 


0.032 


a, cm 


100 


50 


50 


87 


200 




1.6 X 10^ 


4.5 X 10=^ 


7.2 X 10^ 


1.3 X 103 


1.3 X 103 




2.0 X 10^ 


1.6 X 10^ 


1.5 X 10^ 


1.4 X 10^ 


5.5 X 10^ 


Vli , s 


1.7 X 10^ 


9.6 X 10^ 


1.7 X 10^ 


2.6 X 10*^ 


2.4 X 10^ 


te, s 


0.5 


0.1 


0.3 


0.3 


3.5 


e = Pi/a 


5 X 10^3 


1.21 X 10-2 


6.9 X 10-3 


5.2 X 10-3 


2.28 X 10-3 


A 


2.5% 


3.5% 


1-4% 


4.5% 


4% 



Table 1. Typical length and time scales in selected fusion devices (approximate). JET 
parameters estimated from [15], DIII-D from [ ], TFTR from [ ], and ITER from 
TRANSP studies [IS, j'j]. 



Sections 2-10 present a pedagogical derivation all the way from the Vlasov-Landau- 
Maxwell system of equations to the non-equilibrium mean-fleld thermodynamics of the 
system. The order of presentation is the order of the asymptotic expansion, in the spirit 
of asymptotology [20]. However, there are three conceptual strands interwoven in the 
derivation that deserve to be highlighted separately. 

Firstly, there is the equilibrium - the instantaneous solution for the mean flelds. 
First we recover the toroidally rotating Maxwellian for the mean distribution function 
(Section 6.1) and flnd that the toroidal rotation is a rigid-body motion of nested flux 
surfaces (Section 4.1 and Section 6.1). This solution has an arbitrary density and 
temperature for each species and an arbitrary angular velocity on each flux surface. 
Then, in Section 6.3, we flnd that the the poloidal variation of the density is determined 
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by the balance between centrifugal forces and electrostatic fields set up within a flux 
surface. With these results in hand, the poloidal flux function is determined, in 
Section 7.2, by a generalized Grad-Shafranov equation (136). Finally, the first-order 
correction to the mean distribution function is given by the solution of the neoclassical 
drift-kinetic equation, which is also derived within our unified formalism (Section 7.1). 

Secondly, there are the fluctuations that feed on the free-energy sources (gradients) 
present in this (local) equilibrium. The fluctuating distribution function splits into 
the Boltzmann response to the fluctuating electromagnetic fields and a distribution of 
charged rings (Section 6.4). This non-Boltzmann part of the distribution function is 
governed by the gyrokinetic equation (Section 7.4). This system of equations is closed 
by Maxwell's equations for the fluctuating electromagnetic fields in Section 7.5. Finally, 
in Section 9.2, we conclude this strand by deriving the conservation law that governs 
the nature of the fluctuations: the conservation of free energy. 

Thirdly, there is the long-time evolution of the mean fields. This starts in 
Section 7.3, where we determine the evolution of the mean magnetic field. This turns 
out to be completely independent of the fluctuations. The back-reaction of the small- 
scale turbulence on the mean profiles is found in the transport equations of Section 8. 
Examining particle, momentum and energy conservation, we find (166), (179), and (194), 
which determine the transport-timescale evolution of the density, angular velocity and 
temperature profiles in terms of fluxes and sources, which in turn are given as functions 
of the turbulent and neoclassical distribution functions and fields. This strand concludes 
with the results of Section 9.1 and Section 10: that the total energy is conserved on the 
transport timescale and that the increase of mean entropy can be written in the usual 
way as a combination of heating terms and the product of fluxes and thermodynamic 
gradients. 

This paper is written in an entirely self-contained way and so presents both 
rederivations of many known results, cast in forms suitable for our unified framework, 
and a number of new results - we have, throughout the exposition, striven to give credit 
where credit is due without attempting to provide a fully exhaustive literature review. 
To guide a reader aiming to learn gyrokinetic theory from this paper, it is perhaps useful 
to put our approach into the context of other, alternative, approaches. Our exposition 
falls very much within what might be termed "traditional gyrokinetics," in which the 
theory is viewed as an order-by-order asymptotic expansion. We have made no attempt 
to adjust our equations to achieve a Hamiltonian structure or exact energy conservation 
within each asymptotic order (see [21] and references therein) - indeed, the salient point 
of Section 9 is that energy in a multiscale system flows between fluctuations and mean 
fields and so "between different orders" of the asymptotic expansion. The expansion 
terminates at the transport order (Section 8) in the sense that the equations are closed 
and energy is conserved overall. Similarly, we make no attempt to formulate "global" 
equations that simultaneously describe the long and short scales - in fact, we take scale 
separation to be a virtue and consistently enforce it within our formalism. This said, 
alternative approaches (sometimes termed "modern gyrokinetics" ) have many virtues of 
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their own to recommend them - not least some fascinating mathematics - and we refer 
the curious reader to recent reviews [21, 22] where they are presented. 

In Section 11, we present the low-Mach-number hmit of the system of equations 
given in this paper. This hmit removes many cumbersome technical complications and 
so Sections 11.2-11.6 can be used as a concise summary of the basic structure of our 
multiscale hierarchy. It is this set of equations which is implemented in current linked- 
fiux-tube transport codes [23, 24]. 

Finally, in Section 12, we finish the paper by summarizing the conclusions of this 
work and how they fit into the broader landscape of fusion plasma physics. 



2. Fundamental Equations 

As our starting point we take the Fokker-Planck kinetic equation for fg, the distribution 
function of species s, 

where Zg is the charge of the particles of species s as a multiple of the fundamental charge 
e, rris their mass and v their velocity. We will work in Gaussian units throughout with 
c the speed of light, E the electric field and B the magnetic field (throughout this 
work, tildes denote exact fields containing both the mean and the fiuctuating parts; see 
Section 3.1). On the right hand side, C[f] is the Landau collision operator and Ss an 
arbitrary source term that stands in for all physical processes not yet accounted for, e.g. 
atomic physics, fusion reactions, Bremsstrahlung, radio frequency heating and current 
drive. 

The electric and magnetic fields obey Maxwell's equations: 

V-E =Atiq, (2) 
V-S =0, (3) 

f--cVx^, (4) 

^ ~ An- IdE 

where the charge density 'q and current j are 

Q = Y,Zse f d'vfs, (6) 



j = J2Zse d'vvU (7) 
For the purposes of this work, we neglect Debye-scale effects and relativistic effects 



I.e. 

.2 \2 



kixL « 1, (8) 

"-f « 1> (9) 
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where Aoe = a/Te /Airnee"^ is the Debye length, Vth^ = a/ 2Ts/ rus the thermal speed and 
Ts and are the mean temperature and density of species s, where the notion of mean 
will be rigorously defined in Section 3.1. The assumptions (8) and (9) are well satisfied 
in most modern fusion experiments and will be in future ones. 

We use (8) to replace (2) with the quasineutrality constraint | 

^=0, (10) 

and (9) to drop the displacement current in (5) (the Darwin, or quasi-static, 
approximation), giving Ampere's Law, 
~ An- 

VxB = —j. (11) 

c 

We also satisfy (3) and (4) by introducing the scalar potential ip and the vector potential 
A, 

E^^V^~l§, (12) 

B = V X A. (13) 

We will work in the Coulomb gauge V- A = 0. Thus, the four Maxwell equations (2)-(5) 
are replaced by (10), (11), (12) and (13). 

3. Fluctuations and Mean Fields 

3.1. Small-Scale Averaging 

In order to analyse the mean and fluctuating quantities separately, we introduce the 
concept of an average over fluctuations, denoted by Formally, we demand that it 

separate any arbitrary physical quantity g into an averaged part {g)^.^^,.^ and a fluctuating 
part 6g = g — (fi')turb' which by construction vanishes under the average, {Sg)^^^^.^^ = 0. 

As there is a separation of scales, we can interpret this average as an average over 
the fluctuating temporal and spatial scales. The equilibrium length scale, denoted a (and 
understood to be, e.g., the tokamak minor radius), is well separated from the fluctuation 
length scale, taken to be the Larmor radius ps = fth^/^s, where fls = ZgcB / rrtsC. 
Therefore, we can pick an intermediate scale A that satisfies 

a > A > (14) 

and define the perpendicular spatial average of a function g{r,v,t) by 



J d^r'^g{r'^J,v,t) / J (fr^ 



{g{r,v,t))^= d'r'^g{r'^,l,v,t) / d^r^, (15) 



I The quasineutrality constraint implicitly defines the fluctuating electrostatic potential and the mean 
electric field within a flux surface but cannot be used to determine the mean radial electric field [25], 
so we derive an equation for the radial electric field from momentum conservation in Section 8.2 . 
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where A^^ is a small surface which is everywhere normal to the magnetic field and has 
spatial extent of the order of A in both perpendicular directions. The integrals are taken 
at constant v, t and /, where / is the distance along a given field line, or any other field- 
aligned coordinate - see Fig. 1. Clearly, an averaged function cannot vary on the small 
length scales ~ ps, and any function that varies only on the equilibrium length scale 
~ a is unaffected by the average. Similarly, the typical fluctuation time scale u!~^ and 




Figure 1. A section of a toroidal flux surface showing field lines and, in blue, the 
small perpendicular patches A^^ over which fluctuations are averaged in (f5). Black 
arrows denote the normals to the patches, and are aligned with the magnetic field. 

the timescale of the evolution of the mean profiles, taken to be the transport time te, 
are also well separated. Therefore, we can pick an intermediate time T that satisfies 

> T > uj-^ (16) 

and use it to define the temporal average by 

t+T/2 

{g{r,v,t))^ = ^ j dt'g{r,v,t'), (17) 

t-T/2 

with the integral taken at constant r and v. Once again, we see that averaged functions 



Multiscale Gyrokinetics for Rotating Tokamak Plasmas 



9 



cannot vary on the short timescale oo~^, and that functions that vary on the transport 
timescale are unchanged. 

We now define the full average over fiuctuations by 

{9{r,v,t))^^^^ = {{g)^)^. (18) 

We can split the distribution function fs and all fields into mean and fiuctuating parts: 

E = E + 5E, 
B = B + 5B, 
A = A + SA, 

Lf = Lf + Sip, 

Let us use the average (18) to restate the formal assumptions about temporal and 



spatial variation of the mean and fiuctuating quantities ((5')turb ^d)'- 

^ln5(? -a;, (25) 

V\Yi{g),^^^ ^b-V\n6gr~.a-\ (26) 

V±\n5g r^k^r^p;\ (27) 



where b is the unit vector in the direction of the averaged magnetic field, and _L and || 
denote components of vectors or operators perpendicular and parallel to the averaged 
magnetic field, respectively. In Section 3.4, we will formally order all time scales, spatial 
scales and fiuctuation amplitudes with respect to the small parameter e = ps /a. 



Fs 
E 

B 

A 

f 



\/s)turb' 

e' 



B 



A 



turb 



turb 



turb 



turb' 



(19) 
(20) 

(21) 

(22) 
(23) 



3.2. Axisymmetry and Magnetic Geometry 

We will assume the axisymmetry of all mean quantities. In the case of modern tokamaks, 
the deviation from axisymmetry in the magnetic field is extremely small, much smaller 
than Ps/a in such devices. Let us introduce the cylindrical coordinate system: major 
radius R, vertical position z and toroidal angle (p. We pick this system such that {R, (p, z) 
is a right-handed coordinate system - see Fig. 2. Assuming that the toroidal magnetic 
field variation is the same as the toroidal variation of any averaged quantity, we formally 
demand that for any quantity g{r), 



0. (28) 



At the end of Section 5.2, we will extend the notion of axisymmetry to distribution 
functions, which depend upon velocity space as well as the position r. The averaged 
form of (13) is 
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where 

A = A^V(f) + ArVR + A,Vz (30) 

and we have used (28) to drop all (j) derivatives. Therefore, the magnetic field can be 
written in the usual toroidal decomposition [26]: 

S = JV0 + V^xV0, (31) 

where 

i){R, z) = A^ = R^A ■ V0 (32) 
is the poloidal flux function and 

^c'-' - - ^) ^ ■ '''' 

The toroidal symmetry guarantees the existence of well-deflned flux surfaces [27, 28]. 
Topologically, these are nested tori. Since B ■ Vip = 0, these surfaces can be labelled by 
ip. We will see that many mean quantities will only depend on R and z through ^(-R, z). 



3.3. Flux-Surface Averaging and the Motion of Flux Surfaces 

It will be convenient to deflne an average over the surface labelled by tp, which we do 
as follows. For an arbitrary function g{r), 



(gir))^ W = lim^ <j J d'rg{r) j J d'r } , (34) 

A(V) 

where the domain of integration A (■?/') is the annulus between the flux surface labelled 
by ip and that labelled by + Ax(j (see Fig. 3). 
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Note that 

A f (frg = lim I /" d^rg - f d^rg 



lim [ d^rg = V'(g) 



(35) 



where D[il>,t) is the volume enclosed by the flux surface labelled by ip (i.e., the volume 
between that surface and the magnetic axis; see Fig. 3), V = /□(^t)'^^^ is the volume 
of this region and 

If.. dV 



V = lim / ^ , , 

AV^O All) J A (Jll> 



Thus, an alternate deflnition of the flux-surface average is 

IK 1 -9 /• ^3 Id 



\ f dS 
dip / ——g 

JdD(^',t) Nvl 



1 



dS 



(36) 



(37) 



where the surface integral is taken over the boundary of D, dD{ilj,t), and we have 
defined ip in such a way that = at the magnetic axis.§ 




^|J + A^lJ 



Figure 3. A section of tlie torus sliowing tlie regions D{tp,t) and A{ip,t) used in 
defining the fiux-surface average (34), (37) and the surface dD separating these regions. 

§ In our choice of coordinates, ip increases away from the magnetic axis if the toroidal current flows in 
the negative direction and the toroidal field is in the positive direction. Similarly, if the toroidal 
field and current are both in the positive </> direction, ijj will decrease away from the magnetic axis. If 
the positive (j) direction is fixed to be opposite to that of the toroidal current then ij; will always increase 
outwards. 
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We will discover in the subsequent sections (Sections 7.3, 8 and 9) that flux-surface 
averaging allows us to close evolution equations for the mean fields. Two mathematical 
identities will be useful in those derivations: the formula for the fiux-surface average of 
a divergence and that for the flux-surface average of a time derivative. 

First, using (37), we find for any vector field A(r), 



Therefore, 



(V • A), 



(V ■ A)^ 



1 d 



£rV ■ A 



dS 



V {A ■ Vz^), 



A ■ W- (38) 



(39) 



which is the first of the two identities we will require later. This identity also implies 
that 



{B ■ Vg)^ = (V ■ {gB)),^ 



1 d 



V {gB ■ V^P)^ = 0. 



(40) 



'V \ V'dip 

Thus, the fiux-surface average annhilates the operator B ■ V [29]. 

Since the constant-?/' surfaces change in time, fiux-surface averages do not commute 
with time derivatives. Let us consider 



d_ 
Ft 



(frg 



,3 9g 



dD{i,,t) 



dS 



gV^ ■ 



D(V,t) 



(41) 



where is the velocity with which the boundary oi D{ijj, t) moves. The time derivative 
in the left-hand side of (41) is taken at constant flux label ip. Taking the derivative of 
(41) with respect to ip and using (37), we find 



V'Fp 



>3 9g 



V' dt 



Finally, as the boundary of D{ilj,t) is defined to be a constant-?/' surface, we have 



dt 



+ v^-viJ = o. 



(42) 



(43) 



As is defined to be the velocity of a flux surface, we can demand that has no 
component in the surface and so 

dip Vip 



Therefore, 



di\ 



dt |V?/'|2' 

_ ^ d_ 
rV'dt 



V'Fp 



V'(g 



(44) 



(45) 



which is the second of the identities we were seeking. 
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(46) 



3.4- The Gyrokinetic Ordering 

To proceed further we need to impose order on the three small parameters we have, 
e = Ps/a, oj/Qs and the anisotropy of the turbulence k\i/k±, as well as on the amplitudes 
of the fluctuations. Then each term in the kinetic equation (1) will have a well defined 
order in terms of e when compared to flsfs- 

We postulate the following standard ordering [10] of the time and space scales with 
respect to e[|, 

\6B\ \6E\ 6fs u Ps 

\B\ \E\ fs kj_ ^Ig a 
Note that when treating the electromagnetic fields, we order the electric field E and the 
magnetic field B as 

E ~ (47) 

c 

which is equivalent to assuming that the E x B flows are at most sonic and not 
relativistic. 

From these assumptions, we can make a simple random-walk estimate of the 
turbulent thermal diffusivity XTs ~ (gyro-Bohm diffusion [30]) and then order the 
transport time on the basis of this estimate: 

1 XT w /'ps\'^, 

Te a? V ( 

For the collision operator C[/s], we choose an ordering that allows the plasma to 
be either coUisional or coUisionless, namely C[fs\ ~ oofs- Thus, all collision frequencies 
V are ordered 

k ~ 

and any further assumption about the coUisionality will be handled as a subsidiary 
expansion. Even though the collision frequency u is often smaller than w, the collision 
operator C[5fs\ must be retained as the turbulence will otherwise generate arbitrarily 
small scales in velocity space [11, 13, 14, 31, 32] ^. We will return to this point in 
Section 9.2 when discussing free-energy balance. 

Finally, we order the source term as Sg ~ FsIte which restricts us to sources 
of particles, energy and momentum that do not alter the Maxwellian form of the 
equilibrium distribution function (found in Section 6.1). In [ ], this restriction is 
relaxed to consider sources that produce high-energy tails and other non-Maxwellian 
distributions. 

II In this paper, we use the symbol ^ to mean "is the same order as" rather than the more usual "is 
asymptotically equivalent to" . 

^ We choose not to order velocity-space derivatives, which would allow the introduction of smaller 
collision frequencies whilst retaining C[fs] ^ i-^fs- 



~ e^^],. (48) 
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We can now expand the mean (large-scale) distribution function Fs and the 
fluctuating (small-scale) distribution function 6fs as follows 

F, =Fos + F,, + F2s + ---, (50) 
Sfs = Shs + + ■ ■ ■ , (51) 

where Fqs ~ fs, Fu ~ Sfu ~ efs, F2s ~ 5/2s ~ e^/s, etc. 

In the following sections, we insert (50) and (51) into the Fokker-Planck equation 
(1) and expand order by order in e. 



4. Zeroth Order O {Vtsfs) 

In this section, we derive the lowest-order implications of imposing the ordering of 
Section 3.4 on the equations of Section 2. 

We start by taking the lowest-order components of (10) and (11). The former just 
states that the densities must satisfy the quasineutrality constraint, 

Z.eus = 0, ns= j £vFos. (52) 



Substituting the toroidal decomposition of the magnetic field, (31), in Ampere's law, 
(11), we merely find that the mean current j = (j) must be zero to this order, 

\ / turb 

j = ^ Z.eUsUs = 0, Us = — I d^vvFQs. (53) 



By taking the lowest-order component of (12), we also learn that the mean electric 
field is predominantly electrostatic: 



E 



4-1. Sonic Flows 



.Vy, + 0(e2^5). (54) 



We now prove that if the plasma has any perpendicular flow faster than the drift velocity 
then it is a purely toroidal E x B flow. This restriction on the flow has previously been 
found in the context of neoclassical coUisional transport [34, 35, 36, 37]. 
From (1), we flnd to lowest order in e, 

1 \ OF 

E + -V X B] ■ = 0. (55) 
c J ov 

Multiplying by v and integrating over all velocities, we obtain 

E + -u,x B = 0. (56) 
c 

This implies that the perpendicular part of Ug is species independent. Using (56) and 
(54) we have 

E-B = B-Vip = 0. (57) 
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Therefore, to lowest order, the electrostatic potential is a flux function and we can 
partially solve for ip as follows: 

= $ (^) + (58) 

where ipo ~ e$. There is some arbitrariness in the definition of ipo by (58) as we can 
add any function of to it. We resolve this by requiring that <^o vanish when averaged 
over a flux surface: 

(^o)^ = 0, (^)^ = $ (^). (59) 

Solving (56) for Ug shows that any sonic flow that is present is a E x B flow plus 
some arbitrary parallel flow 

Us = —b X V(p + Uiub, (60) 
13 

where b = B /\B\ . Using the solution (58) for (p and the toroidal decomposition of the 
magnetic fleld, (31), to expand b, we have 



Us = u{ip)R'^V(f) + 



uns-uj{ip)^ 



= (61) 



The flrst term in (61) is a species- independent rigid-body toroidal rotation of each 
individual flux surface with an angular velocity (^{ip). This will be denoted by 

u = uj{i))R^V<p. (62) 

The second term is a purely parallel flow, which in Section 6.1 will be shown to vanish 
to lowest order. 



5. Rotating Gyrokinetic Variables 

Before continuing to expand the kinetic equation (1), it will be convenient to introduce 
the following new variables [1, 10, 35, 38, 9]: the guiding-centre position ii^, particle 
energy e^, magnetic moment /i^, gyrophase and the sign of the parallel velocity a. 
The variable transformation of the phase space is 

(r, v) {Rs, Es, Us, 1^, cr) (63) 

and the new variables are deflned by 

„ b X w , , 

Rs = r- (64) 

Ss =^msV^ + Zse[mj) + ipo]-Zse<^irs), (65) 

P^s = (66) 

a (67) 
Fill 

where w is the peculiar velocity of the particles with respect to the toroidal rotation: 
w = V — u = w\\b ^ w± (cos "i? 62 — sin-i^Ci) , (68) 
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with the toroidal velocity determined by (62), Ci and 62 arbitrary orthogonal unit 
vectors perpendicular to the magnetic field (with b = 62 x ei) + , and ip* is a fiux-like 
quantity proportional to the toroidal canonical angular momentum: 



^PKr^v) = ^(r) +'4^{vV<f>)R'=^ + '^{w- V0) R' + (69) 



The new velocity variables Eg and /i<j are closely related to conserved quantities of 
the particle motion in the mean electromagnetic fields. The magnetic moment Hs is the 
first adiabatic invariant of the particle gyromotion [40, 41, 42]. The energy variable. Eg, 
is constructed from the conserved total energy and the conserved quantity ips so that 
it is the energy in the rotating frame to lowest order. This can be seen by expanding 
^{ipD around yielding to lowest order 

Es = ^rrisW^ - ^msU^{i')R'^ + Z^eipo + 0(erj. (70) 
The kinetic equation (1) can be written in these new variables as 

dfs dfs ■ dfs dfs dfs -dfs ^r.n , ^ 

where g = dg/ dt denotes the time derivative of g along a particle orbit. The expressions 
for Rs, fis, £s, and are derived in Appendix A, see (A. 9), (A. 25), (A. 33), and (A. 39). 



5.1. The Gyroaverage 

Transformation from the spatial coordinate r to the guiding-centre position Rg (known 
as the Catto transformation [')]) allows us to introduce gyroaveraging — a crucial 
mathematical device that will enable us to close the equations for and 6fs- We 
define the gyroaverage of a quantity g by 

{9)r = ^ ^ d'd9{Rs,£s,fJ^s,^,(T), (72) 

where the integral is performed holding Rg, Eg, [is constant. 



5.2. Derivatives and Integrals 

All quantities that have velocity-space dependence will be functions of Rg, Es, Us, 
and a. The quantities that are only functions of space will be evaluated at the spatial 
position r, unless explicitly stated otherwise, with r considered as a function of Rs, [ig, ^. 
In what follows, V is reserved for derivatives with respect to r, while derivatives with 
respect to Rg will be written explicitly as d / dRs . 

As the variables we have chosen are adapted to the particle motion, we will not 
transform the field equations into these variables. The fields are functions of space 
but not velocity so they remain functions of r and t in accordance with the principle 

+ There might be a concern that we cannot define such an ei and 62, this is resolved in [39] where 
it is proved that for toroidal confinement devices we can always make a globally-valid choice of these 
vectors 
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stated above. Therefore, the integrals over velocity in the definitions of density, (6), and 
current, (7), are to be taken at constant r, so functions of Rg, Eg and fig appearing in 
the integrand are to be considered as functions of r and v via the definitions (64), (65) 
and (66). Namely, for any function g, 

jd^wg{Rg,es,Hs,^,(T) = jd^wg{Rg{r,w),eg{r,w), ^g{r,w),i^{r,w),a{r,w)). (73) 

This is how all velocity-space integrals will be performed in this paper unless explicitly 
stated otherwise. With this definition, integrals over v and w at constant r are 
equivalent. 

Finally, as promised in Section 3.2, we now give a precise mathematical formulation 
of the assumption of axisymmetry of mean fields as applied to the mean distribution 
function: 

(^*' -A 

which, to lowest order in e, is equivalent to 

(V0)-V^,^^,,F, = O. (75) 

However, this is not equivalent to (V0) -Vj^F, = 0: indeed in (68), the basis vectors 
b, Ci and 62 possess non-zero variation in the direction. 

5.3. Gyrotropy of Fqs 

To zeroth order, the kinetic equation (71) turns out to be just 

= 0. (76) 

■f^S )/^S I^S 



Fs = 0, (74) 



Thus, Fqs is gyrophase independent. Note that this means that, to lowest order, 

U^ww^Fos = 0, (77) 



confirming the earlier result that the lowest-order perpendicular velocity is completely 
contained in Ug as defined by (61). 

6. First Order 0(eJ7s/s): The Maxwell-Boltzmann Equilibrium 

In this section, we start from the kinetic equation (71), and expand it to first order. 
Using the lowest-order expressions for Rg, ig and t? given by (A. 10), (A. 34), and (A. 39) 
respectively we find 

(u.||b + . ^ + - Zgew^ ■ VVtt^ + fis^ iF,g + 6fu) = C[F,g] , (78) 

^ ' oRg o^ig osg ov 

where 6ip' is the fiuctuating electrostatic potential in the toroidally-rotating frame, given 
by 

S^' = S^--u- 6A, (79) 

c 
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and u without the species index is the toroidal-rotation component of Us defined by 
(62). 

In the following subsections, we analyse (78) to discover that: the lowest-order 
mean distribution function Fos is Maxwellian; the bulk motion is the purely toroidal 
rotation of flux surfaces, i.e. the second term in (61) vanishes and Us = u; and Sfg 
consists of the Boltzmann response to 6ip' and a gyrophase-independent distribution of 
charged rings hs{Rs, /^s, e^, cr, t). 

6.1. Maxwellian Equilibrium 

We flrst prove that Fqs is Maxwellian. Gyroaverging (78) and using the fact that, to 
lowest order, {^s) r = (see (A. 29)), {w±_ ■ W6(p')j^ = (see (A. 4)), and u-dFog /dRg = 
(axisymmetry), we obtain 

^\\b-g£ = {C[Fos])R. (80) 

Multiplying this by 1 +lnFos, integrating over all velocities and averaging over the flux 
surface, we flnd 

j d'ww^i b ■ ^ (Fo, In Fo,) ^ =(^j £w In Fq, C [Fq,] 

To lowest order, we can replace all instances of Rs on the left-hand side of this equation 
with r and write the velocity-space integral in terms of integrals over Eg, fJ^s, and a at 
constant r using 

J mi\w\\\ 



This gives 



2^E^/ ^^■V(^o.lnFo.)\ (83) 



V ■ ( b ld^ww\\FnMF< 



By using (39) to express the action of the flux-surface average on a divergence and the 
fact that b ■ Vip = 0, we conclude that the above expression vanishes. Therefore, from 
(81), 

jd'w\nFosC[Fog]^ =0. (84) 
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By Boltzmann's if-Theorem, (84) implies that F^g is a local Maxwellian [43]. We know 
that, by definition, this Maxwellian has density Ug, temperature Tg and velocity Ug, 
where Ug is given by (61).* Thus, Fog can be written as 



Ugir 



rrig 



3/2 



nig 



exp 



- 2mgW\\u\\g{r) + u\g{r) 



2Tg{r) 



where 



till 



Mil 



B ■ 



^5) 



(86) 



This form does not contradict the condition (76) that F^g must be gyrophase- 
independent whilst holding Rg, Sg and Hg fixed because the difference between Rg and 
r is higher order. However, we still wish to have F^g expressed in the [Rg, Eg, fig, a) 
variables. We accomplish this by using (70) to express irtgw'^ in terms of Eg and find 

Eg mgW\\U\\g{R, 



Os 



Ng Rg 



rrig 


3/2 


exp 


_27rTg{Rg)_ 



Tg (Rg 



Tg(Rg 



+ 0{eFg), 



7) 



where 



g = Ug exp 



mgU:'^{il))R^ ZgCifo 



2Tg 



+ 



Tg 



+ 



2Tg 



and we have used the fact that to lowest order the mean fieds taken at the guiding-centre 
position Rg are the same as when taken at the particle position r. 

Inserting F^g, given by (85), back into (80) and dividing through by F^g we get 

3 



w\\b-V \ \\iN. 



InT, 



W\\Eg 

-j^b-VTg-w\\h-V 



mgW\\u\\g 



0. 



9) 



This equation must hold for all velocities w, so each term in this equation must vanish 
independently. Tackling the second term first, we see that the temperature Tg must be 
a fiux function to lowest order: 



b-VTg = 0, 

so Tg = Tg{ip). For the first term to vanish the same must be the case for Ng 
b ■ VNg = 0, 

so Ng = Ns{iIj). Turning now to the third term of (89), we see that 



2wib ■ Villi 



u 



\gb ■ Vwy 



0, 



(90) 
(91) 
(92) 



which can only be solved for abritrary w\\ by u\\g = 0. Thus, the background Maxwellian 
only depends on Eg and so is isotropic in w and the-lowest order (sonic) fiow is a pure 
toroidal rotation: Ug = u = uj{il))R^'W (j). As this fiow is species-independent there are 

* Formally, (84) implies that all species have the same temperature and mean velocity. We will show 
that they do indeed have the same mean velocity, but we will gratuitously retain the species dependence 
of the temperature. This will only be of importance when we come to discuss heat transport and so 
we defer the discussion of interspecies temperature differences to the end of Section 8.3. 
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no currents in F^g, since the plasma is quasineutral (52). This is consistent with the 
lowest-order Ampere's Law (53). The vanishing of the parallel component of (61) is in 
accord with the well known result [36, 43, 44] that poloidal flow is strongly damped 
on the ion-ion collision time, which, in our ordering is, indeed, 1/e^ shorter than the 
timescale of the evolution of the bulk flow. 

We can finally write the complete solution for Fqs as 

-1 3/2 



-es/Ts{iP(Rs)) 



(93) 



_27rT,{^{Rs)) 

This form of the distribution function is manifestly gyrophase independent holding R 
and Eg constant and so is consistent with (76). 



6.2. Gyrotropy of Fig 

Substituting (93) back into the first-order kinetic equation (78), we find 

^s— (Fi, + 5/i,) = -^w^ ■ Vd^'Fos. 
ov Ts 

Averaging this over the fluctuations gives 
dFu 



0. 



(94) 



(95) 



Thus, Fis is gyrophase independent. 



6.3. Poloidal Density Variation 



Rearranging 
given by 



and using u\\s = 0, we see that the physical density of a species is 



Hs = Nsiij) exp 



+ 0(ens 



(96) 



2T, T, 

which is not a flux function, unless the rotation is subsonic (see Section 11) [36, 43]. 
We can insert this into the lowest-order quasineutrality condition (52) to determine ipo, 
subject to the constraint that (v^o)^ = (by definition, see (59)): 



ZsNs{ip)exp 



2T, 



0. 



(97) 



Physically, this demonstrates that the poloidal variation of the density on a flux surface 
is governed by the balance between centrifugal forces which sweep heavy particles (i.e. 
ions) to the outboard side (higher R) and the electrostatic potential that is set up within 
the flux surface to maintain quasineutrality. 
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6.4- The Boltzmann Response 

Taking the lowest-order fluctuating components of (12) and (13), we see that the 
fluctuating electric field in the toroidally rotating frame is electrostatic, 

6E + -ux6B = -V±6<^' + 0(e^E), (98) 
c 

where 6ip' is defined in (79). Using (95), the first order kinetic equation (94) is 

^s^^ = -^w^.V6^'Fos. (99) 

OV Is 

We can integrate this with respect to i), using (A. 3) for 6(p': 

d6if' 



and obtain 



(100) 



6fu = -^^^^^Fo, + hs {Rs, /i„ Bs, (T, t) , (101) 

s 

where the constant of gyrophase integration hg is the gyrophase-independent distribution 
of Larmor rings, which will completely describe the kinetics of the small scale turbulence. 
The gyrophase- dependent part of Sfis is just the Boltzmann response in the moving 
frame to the electrostatic field (98). 

6.5. Summary: The Lowest- Order Solution 

Collating the results of the previous subsections we have the solution for fg to first order 
in e: 

fs =Fs + 6fs, (102) 
Fg =Fog{ij{Rg),es) + Fu{Rs,es,fis,(T) + 0{e'f), (103) 

^fs = - -TfSip'{r)Fos + hs {Rs, Bs, fis, a) + (eV) , (104) 



Fos = Nsii^iRs)) 



3/2 



g-../T4V(fl.))^ (105) 



_27TTs{ij{Rs))\ 

In the next section we find equations that govern the gyrophase-independent functions 
Fls and hs- 

7. Second Order O (e'^^lsfs)'- Neoclassical Theory and Gyrokinetics 

In order to completely determine the distribution function fs given by (102), we need 
equations for Fis and hs. In order to find them, let us substitute the form of fs 
summarized in Section 6.5 into the exact kinetic equation (71) and keep only terms 
up to second order in e: 

-^SQ^ {F2S + Sf,s) + ^ ( -^^Fos ] + C[Fos + Fls + hs] , 
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where d/dt is the full derivative along a particle orbit as in (71) and we have used the 
gyrophase independence of Fq^ (76), Fis (95) and hg (101). Equation (106) contains 
and Sf2s, about which we have no information. It is thus crucial that Fqs, Fis and hg 
are gyrophase-independent, as this allows us to gyroaverage (106) to find an equation 
that does not contain any second-order distribution functions: 



dhs 
dt 



+ (R 



d 



R dRs 

Fo, I d 



R 



+ 



dt 



{Fos + -^Is 



hg 



+ {C[Fos + Fu + hg]) 



(107) 



R 



where we have used Fqs given by (105) and the fact, derived in Appendix A, that 
{[is)r = Oie^nsfis) (see (A.29)) and {is)^ = 0(e2^],r,) (see (A.35)). 

The gyroaverages appearing in (107) are calculated in Appendix A: we use (A. 23) 
for (^Fts^ and (A. 55) for the combination of gyroaverages appearing on the right-hand 
side: 

d 



d 



+ u(Rs 



d 



dt 

ZsGFqs 



d_ 
dt 

11 



+ uiR, 



■ 

d 



d 



OR. 



dRs 

{x)r 



(Fi, + hs) + {Vus + ■ ^ (Fos + hs) 



rrisFps 
Ts 



Iw\\ 



+ to{ij)R' 



du 
dip 



'Vx)r-^^ (108) 



Z e 3 A 

■'-w,Fos—-b+{C[Fos + Fu + hs]) 



dt 



where we have defined the gyrokinetic potential 



1 



c 



V ■ SA = b^p w ■ SA, 

c 



the associated fiuctuating velocity field^ 

c 



B 



b X Vx, 



R 



B 



b X 



d{x)n 
dRs 



(109) 



;iio) 



and the guiding-centre drift velocity 

b 



Vy,s 



ns 



wtb-Vb+ -wiVlnS 



-uj\il))RVR - 2w\\uj{il))b X V2 + — V^o 



nrts 



:iiii 



which consists of, in order of appearance in (111), the curvature drift, the V-B drift, the 
centrifugal drift, the Coriolis drift, and the mean first-order E x B drift. 

The detailed derivation of (108) is given in Appendix A. 6. In the rest of this section, 
we will use (108) to derive a closed solution for Fis and hs, accurate to first order in 
e, in terms of the unknown functions Ns{ip), T,, u}{ip), I and ip, which parametrise the 
equilibrium. Equation (108) can be split into a mean equation (see (112)), which will 
determine Fis, and a fiuctuating part (see (144)), which will determine hg- We will find 

tt This can be shown to consist, physically, of the fluctuating E x B drift in the rotating frame, the 
motion of guiding centres along fluctuating field lines and the fluctuating VB drift. This is proved in 
detail for gyrokinetics in a non- rotating slab in [ '"]. 
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that the mean equation is precisely the neoclassical drift-kinetic equation [3n] and the 
fluctuating equation is the gyrokinetic equation for a rotating plasma [1, 16]. These two 
equations are closed by the mean and fluctuating Maxwell equations, which relate ijj 
and I to Fis (Sections 7.2 and 7.3) and the fluctuating fields to hg (Section 7.5). 



7.1. Neoclassical Distribution Function 



Averaging (108) over the fluctuations, we obtain 



'C[F,.,]) 



R 



-Vr 



Ds 



OR, 



w\\Fos — ■ b 



dt 



\c[F,;\) 



R- 



;ii2) 



This is just the usual neoclassical drift-kinetic equation [35, 43, 36]. The coUisional term 
involving Fis can be split from that on Fqs because the collision operator is linearised 
around the Maxwellian part of Fq^. Note that C[Fos] is only zero to lowest order, so in 
(112) it represents collisions on the first-order departures of Fqs from a pure Maxwellian 
(having to do with the absorption into Fq^ of some non-Maxwellian velocity dependence 
via ip{Rs) and e^; see Section 6.1 and Appendix D). 

We now wish to solve (112) by first separating particular solutions corresponding 
to some of the source terms in its right-hand side. To deal with the first of these terms, 
we let 



where 



K = 
w\\b ■ 



Fis{Rs, £s, f^s,o-) + F*^ 



nisC 



dFl 
OR., 



B 



dFos 



d 



dR.. 



Iw\\+uj{tl))R^B' 



(113) 
(114) 



dF< 



Os 



'Ds 



v^) 



dF< 



Os 



dip 



-Vds 



dip 
dFos 



;ii5) 



dRs 



where we have used the well-known identity [36] 

7u7|j BR^ 



Vds ■ VV" = -w\\b -V 



+ 



fis 



uj{iP) 



Therefore, 

"^11^ -m: 



c 



Fu 



R 



Z p d A 

~w\\Fos—-b+{C[Fos + Fl2)^. 



;ii6) 



;ii7) 



We now wish to replace the term depending on A with a source term that depends 
on {E ■ B)^, which will later turn out to be opportune (see Sections 7.2 and 7.3). To 
do this, we follow [1, 5, 47] and introduce the ansatz 



Fi, = Fi-)(il„ /i„ a) - ^Fo, 



dl' B 



{E-B), , idA 



+ 



c dt 



;ii8) 
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where the integral is performed along a fixed field line, and /' is the distance along that 
field line. Note that the second term in (118) carries no current. Inserting (118) into 
(117), we find 



w\\b 



dK 



(nc) 



dRs 
Finally, we split 



w\\c {E ■ B) 



V. 



where Fg^^ and Fs^^^ satisfy 
dFP 



Fo, + {C[Fo, + F*^]) 



R- 



(119) 



(120) 



c[fP] 

_p(nc) 



2\ls -Fos, 

{C[Fos + F*,])^. 



:i2i] 



:i22) 



This pair of equations can now be solved without knowing {E ■ B)^ and then Fs^^^ is 
given as a function of {E ■ B)^ by (120). This will be used in solving for the evolution 
of the magnetic field in Section 7.2 and Section 7.3. 

7.2. Ampere's Law and the Magnetic Equilibrium 

The average of Ampere's Law (11) over the fiuctuations is 

j = ^VxB, j = Yl I d^^wPs, (123) 

where j is first-order (to lowest order, jf = 0; see (53)). Using the axisymmetric form 
for B (31), we can express V x S in terms of derivatives of / and tp: 

j = |^[V/x V0-(A» V0] 

where A* is the Grad-Shafranov operator 

92 



(124) 



1 d 92 

+ 



(125) 



dR^ RdR dz"^^ 

We can now use the complete solution for F^ given by (103) to determine j . This 
is done in Appendix B.l (see (B.17)): 



j=cR \^n,\T.— + 



dlnTg 

dip 



V0 



+ cR^ (^m,nsR^^ co{ij)^V(P + K{ij)B, 



(126) 



and we have defined 



(127) 
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with pj^^^ given by the neoclassical drift-kinetic equation (119). The fact that i^' is a 
flux function follows directly from (119), and is derived in Appendix B.l (see (B.15)). 
Equation (126) is consistent with this, as can be seen by taking the divergence of 
both sides. Note that the perpendicular part of (126) is a statement of force balance 
for the mean quantities. Indeed, taking the cross product of it with B and using 
i?2v0 X B = W, we get 



-JxB = Y^ 



dip 



ZsCipo - -msu'^{tp)R^ + Ps 



dlnPs 

dip 



nXsUsR^ 



(128) 



where we have defined the pressure Ps = UgPg and used (96) for Ns{ip) and (62) for u. 
This is just the balance between the Lorentz force, pressure and inertia. 

We will now take three projections of (124): onto Vip, onto Vip x V0, and onto 
V0. Taking the projection onto Vip first, we have, from (124), 



j -Vip = — (V/ X V0) ■ 

4:71 

and from (126) 

j-ViP = 0, 

whence 

(W X V0) ■ V/ = ■ V/ = 0. 

Therefore, / = I{ip) is a flux function. 

Taking the projection onto Vip x V0, we have from (124) 



j ■ {ViP X V0) 
and from (126) 



dl 



47r/?2 dip 



\ViP\' 



j ■ {Vip X V0) = K{iP)B ■ {Vip X V0) = K{iP)\Vip\'^R-^. 
This gives us the following equation for I{ip) 

Finally, the toroidal component of (124) is just 



47r/?2 



A*ip. 



(129) 
(130) 
(131) 

(132) 
(133) 
(134) 
(135) 



Using the toroidal current from (126) and expressing K{ip) via (134), we obtain 



A*ip = - 4nR^ Us I Ps , , + 



47ri?2 




Zsd^f) - ^insUJ^{ip)R^ + Ps 



d\nPs 

dip 



(136) 
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which allows us to determine ijj as a function of R and z if we know Ng, Ts, cu and /. 
This is the generalization of the Grad-Shafranov equation for a rotating plasma. 

7. 3. Evolution of the Magnetic Field 

The results of the previous section complete the solution, including time dependence, 
for the magnetic field. We can solve (119) for F^"""^ to find K{^) in terms of {E ■ 
and then use (134) to find {E ■ B)^ in terms of I^^j). We can then use (33) to find I{jp) 
in terms of A and so our expression for 

is an evolution equation for A (and thus B). Indeed, taking the time derivative of (33) 
divided by i?^, we find 

Flux-surface averaging and using (39) and (45) we have 
d 



dt 



^V'm{R-% = c^V'{E.B)^, (139) 



where we have used (31) to rewrite the right-hand side in terms of B. 

In the study of tokamak plasmas, it is conventional to work with the safety factor 
qlip) rather than The safety factor is defined in terms of the toroidal flux and 

the poloidal flux ip by 

= = ^ / ^'^^ ■ (140) 

2vr # 27r J^^^^^t) 

where the integration is carried out over the volume interior to the flux surface labelled 
by ip (see Section 3.3). Geometrically, q{ip) is the number of times a field line on a given 
flux surface toroidally winds round the vertical symmetry axis for each poloidal transit 
around the magnetic axis. From (140), the definition of the flux-surface average (37), 
and the axisymmetric form of the magnetic field (31), we can express \E' as 

where V is defined by (36). Therefore, 

^W = ir^V'm{R-'),- (142) 
From (139) and (142) we immediatly find 



di 



Thus, q{ip) evolves via (143), /(V') is determined by (142), and ip{R,z) is then 
instantaneously determined from the generalized Grad-Shafranov equation (136). This 
completes our solution for the evolution of the magnetic field. 
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It turns out that q{tp) can be shown to evolve on a timescale that is even slower 
than the transport timescale, the resistive timescale, i.e., c{E ■ B)^ turns out to be of 
the order of Uthi(?7ie/mj)^/^e^-B^. This result is proved in [18], where we also relate 
the evolution of q{ip) to the conservation of magnetic flux. Indeed, the geometric 
interpretation of q{il)) is consistent with the notion that this longer timescale is the 
timescale on which the topology of the magnetic field changes, even if other properties 
of the field may change more rapidly. 



7.^. The Gyrokinetic Equation 



Turning to the fluctuating component of the second-order kinetic equation (10^ 
findtt 



we 



d_ 
di 



+ u(Rs 



d 



Qs 



d 



OR.. 



R 



u(R., 



d 



dFos 



+ 



rrisFos 



dRs_ 
B 



{x)r 



(144) 



with and Vbs given by (110) and (111), respectively. Explicit forms for the collision 
term {C[hs])j^ are proposed in [49, 31, 13] for various model collision operators, however, 
only the properties common to all such operators will be needed in this work so we leave 
this term in its general form. Finally, we reiterate that this equation is written in terms 
of the variables Rg, Es, fis and t?; thus w\\ is given by 



m, 2 



(145) 



It is the gyrokinetic equation (144) which is solved by the gyrokinetic turbulence 
simulation codes [50, 51, 52, 53] and has formed the basis of a large body of theoretical 
analysis of tokamak micro instabilites and turbulence. | The turbulence is driven by the 

ft This is in agreement with the nonhnear gyrokinetic equation of [ ], but not with that of [ ]. The 
discrepancy arises from a subtlety in the gyroaveraging of u ■ Vx that is explained in detail in Appendix 
A. 6 (see (A. 54) and Appendix A. 6. 2.) 

I In the gyrokinetic literature, there is much discussion of the so-called "parallel nonlinearity" [54, 55], 
which is proportional to 5E^^ {^85 fs/ dv\\) . Some derivations of gyrokinetics include this term in the 
gyrokinetic equation, claiming that it has beneficial numerical properties and an important role in 
turbulent heating and conservation of energy. The gyrokinetic equation (144) derived here does not 
contain the parallel nonlinearity because the latter is higher order than the terms we retain and only 
appears in the derivation of the transport equations. We will (in Section 9) derive conservation laws 
that prove that our equations nevertheless conserve, in the mean, energy and satisfy Boltzmann's H- 
theorem. It is true that in our formalism there is no formal energy conservation law for gyrokinetics on 
the fluctuating timescale (in contrast to versions of gyrokinetics derived from Hamiltonians; see, e.g., 
[56]). This is because energy is exchanged between the fluctuations and the mean fields. We will show 
in Section 9.2 that, on fluctuating timescales, our gyrokinetic equation conserves the free energy^ which 
is the relevant conserved quantity for kinetic turbulence [11, 12, 13, 14] and, moreover, the conservation 
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density and temperature gradients in the equilibrium (through the dF^s/ dip term on 
the right-hand side of (144)) and by the velocity gradient (through the du/ dtp term). 
It is important to note that this equation is only coupled to the neoclassical distribution 
function Fs^^^ through the slow evolution of Ns, Ts, u and /. 



7.5. The Fluctuating Maxwell's Equations 

To solve (144), we require a way of obtaining {x)r from hg. This is provided by the 
fluctuating parts of Maxwell's equations. 

The fluctuating component of the quasineutrality condition (10) is 

s * s 

where 6ip' is given by (79). The integral over velocities is performed at constant r as 
discussed in Section 5.2. We have now made this explicit by introducing the gyroaverage 
at constant r, w\\, and w±: 

{hs)r = d^h,{R,{r,wii,w±,^),e,{r,wi\,w±,^),fi,{r,w±),a). (147) 

The fluctuating component of Ampere's Law (11) is 

V xSB = -V^5A = — ^ [ £w{whs)^. (148) 

The parallel component of this equation gives us the following equation for 6A\\ : 

4:71 



Vi 6 An 



J^Zse jd^ww\\{K)^. (149) 



It turns out to be convenient for various calculations to work with 5B\\ rather than 
5 Ax. Thus, we take the divergence of b x (148) to flnd 

V ■ [b X (V X 5B)] = Vi55|, = —J2 [d^wV± ■ {b x w^hs)^. (150) 

We now use 

dw 

bxw = - — , (151) 

integrate by parts inside the gyroaverage, and use (see (A. 3)) 

dh 



w± ■ Vhs = —^Is 



to obtain 

. 8B\\B 



+ 0{tn,hs) (152) 



: J2 J d'winisW^w^K), = 0. (153) 



of which is crucial for the mean energy to be conserved on the transport timescale. Note that numerical 
evidence ['i4] has shown that the inclusion of the parallel nonlinearity in the gyrokinetic equation does 
not affect the solutions, which is as expected for a term that is smaller than all the others in the 
equation. 
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Once again (cf. (128)), perpendicular Ampere's law has become the equation for 
perpendicular force balance. This eliminates fluctuations such as the compressional 
Alfven wave (see the appendix of [57]), which are not force-balanced. 

This completes the description of the fluctuations, which are now determined by a 
closed set of equations: (104), (144), (146), (149) and (153). 



8. Third Order O(e^^lsfs)'- Transport Equations 



So far we have derived equations for the fast evolution of the fluctuations, instantaneous 
relations between equilibrium quantities, and an equation for the slow evolution of the 
magnetic field. These hold for a given set of three flux functions: A^^, u and Tg. In 
this section, we derive equations for the slow (transport-timescale) evolution of these 
functions. We refer to these as the transport equations. 

To derive such equations, we will need to go to the next order in our expansion in 
e, namely 0{e^Qsfs)- As we wish to maintain the physical interpretation of evolution 
equations for Ug, u and Tg in terms of the transport of particles, momentum, and heat, 
we return to the original r and w variables to begin our derivation. The averaged form 
of (1), written in these variables, is 



dt 



+ (u + w)- VFs + 



du 

as - — ~ [u + w) ■ Vu 



dw 



Sag 



dw 



turb (154) 

= C[Fg] + Sg, 

where we have naturally split the particle acceleration into its mean part a^ and its 
fluctuating part Sag. The mean acceleration is 



a. 



ZgC 
mg 



E 



w 



X B 



rrig 



IdA 

c~dt 



(155) 



where we have used the results of Section 4.1. The fluctuating acceleration is 



6a^ 



ZgC 
nig 
ZgC 



6E + - (u + w) X 5B 

c 



rUg 
ZgC 

TUg 



^ 1 1 
Vx+ -w- V6A + - 

c c 



6ip' w ■ 6A 

c 



d_ 
dt 
1 



+ u-V\5A 



(156) 



w±_ ■ V5A 



where x is defined by (109) and 6(p' by (79). Equation (154) will correctly describe the 
transport-timescale evolution of Fg if we keep all terms up to order 0{e^Qgfg). This 
clearly requires knowledge of the distribution function including 0{e^fg) corrections. 
From the previous two sections, the particle distribution function is 

fs = Fog{ilj{Rg),eg) + Fig{Rg, e,, /i^, a) + F2g{r, v) 



ZgC 

~t7 



6ip'{r)Fog + hg{Rg,es,fig,a) + 6f2g{r,v) + 



(157) 



where the equilibrium distribution function Fqs is given by (93), Fig and hg are obtained 
from the neoclassical drift-kinetic equation (112), and the gyrokinetic equation (144), 
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respectively, and we have absorbed all (as yet unknown) higher-order terms into and 
6f2s- The exact electric and magnetic fields are 

1 dA 1 85 A 

c~dt 

and 



E = E + 5E = -V$ - V<y9o - V(5v9 



c dt 



(158) 



B = B 



5B = I{^ilj)W<p + Wip xW(t) + W X 5A. (159) 

In the expressions for the fields, and ip{R., z) are determined as explained in Sections 
7.2 and 7.3, (p^ is obtained from the quasineutrality condition (97), and 5ip' and 5 A are 
given by the fluctuating Maxwell's equations (146), (149) and (153). 

The solution (157) for fs will need to be completed with information about and 
^f2s- We can obtain the gyrophase-dependent part of to second order by expanding 
(154) to second order in e: 



- ( 6a 



w) ■ VF, + 

d6f^ 



' V(/?o + -— I + (*i + if^) ■ 



m 



dw 



c dt J 
+ C[F,] + 0(e='nj, 



m 

dw 



(160) 



turb 



where we have used 

{w X b) 



dw 



d_ 



(161) 



Similarly, the fluctuating part of (1) can be expanded to obtain the gyrophase 
dependence of 5fs to second order: 



d_ 



5fs 



d_ 
dt 



+ u - V ] 5fs - w\\b ■ VSfs - Stts 



m 

dw 



+ 



/ ldA\ 



m 



w) ■ Vu 



dSfs 
dw 



- 6a. 



d6fs 
dw 



(162) 



By examining the order of the terms on the right-hand side of (160) and (162), we see 
that we only need the particle distribution function correct to 0{efs) in order to evaluate 
the left-hand side correctly to order 0{fls^'^fs)- 

In the subsequent sections, we will take moments of (154) and show that (160) and 
(162) contain sufficient information about the second-order parts of the distribution 
function to close the resulting transport equations in terms of known quantities - i.e., 
the gyrophase- independent parts of and 6f2s do not affect transport. 



8.1. Particle Transport 

To derive the particle transport equation, we integrate (154) over all velocities to find 
dn,. 



dt 



+ V ■ (risUs 



q{n) 



(163) 
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where 

= Id^wSs, (164) 



risU, = jd^wwFs. (165) 

Equation (163) clearly describes all particle transport: rig is created or destroyed by the 
source S'i"'' and advected by the mean flux UgUs. As the first and last terms in this 
equation are 0{e^Qsns), we conclude that we require risUs correct to order ©(e^n^fthj 
to evaluate the divergence. It is clear that the parallel part of the flux will depend on 
the gyrophase-independent part of F2s, which we cannot easily find. In contrast, the 
perpendicular flux will only depend on the gyrophase-dependent piece of ^2^, which we 
can find via (160). 

To reduce (163) to an equation only containing the perpendicular fiux, we apply 
the fiux-surface average, defined by (34). This results in 



1. ^ 
V' di 



1 d 



^/'(^s), + y^g^V'{T.)^ = {Si-%, (166) 

where we have used (39) and (45) to simplify the fiux-surface averages of the divergence 
and the time derivative. We have combined the radial particle fiux§ and the terms due 
to the motion of the fiux surface into 

Ts = nsUs-V^lj + ns^. (167) 

We can write the first term of (167) as 



n 



,Us-Vip = jd^w {w^-Vip)Fs = jd^wR^B{v - Vet)) 



(168) 

r,wn ,uix 



where v is just shorthand ioi u + w and we have used 

Fin 

■ V(j) (169) 



w^-Vij = -F?Bw ■ (b X V0) = R^B (bxw)-V^ = -R^B — 

ov 

r,iii|| ,wj_ 

and integrated by parts with respect to {}. So, in order to calculate the radial fiux to 
second order, it is sufficient to know dFg/ dd up to second order only. 

In Appendix C.l, we perform the explicit evaluation of {usUs ■ V'?/')^ via the kinetic 
equation (160), resulting in (C.IO), which we substitute into (167) to find: 



[d'w{h,V^)^-vA ) . 

J I turb/ ^ 



(170) 



The first term in the above expression is due to classical coUisional transport with Fqs 
given by (93) || , the second and third terms are due to neoclassical transport with Fs^'^^ 

§ Whilst this flux is not in the geometrically radial direction, but is in fact the cross-flux-surface flux, 

it is both convenient and conventional to refer to it as the radial flux. 

II Note that F^g is only Maxwellian to lowest order, so C[i^os] = ^{^^^sPqs) is non-zero. 
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and {E ■ B)^ calculated as explained in Sections 7.1-7.3^, and the final term gives the 
turbulent contribution to the particle flux in terms of hg, which is the solution of the 
gyrokinetic equation (144) + . We note that the composition of a fiux-surface average 
(defined by (34)) with the fluctuation average (deflned by (18)) is 

t+T/2 

(171) 

t-T/2 Aa 

where is the annulus bounded by two flux surfaces a distance A apart centered 
on ipif)- This is exactly the average over a flux tube that is implemented in current 
gyrokinetic linked-flux-tube codes [23, 24]. 

Thus, we now have an equation, (166), that determines (^s)^, but we still need to 
use this to flnd Ng and hence via (96). This is done by flux-surface averaging (96) to 
flnd 



exp 



2T, 



(172) 



which allows us to express A^^ in terms of (n<j)^ and known quantities. In particular, we 
can solve (172) simultaneously with the quasineutrality condition (97) to flnd A^^ and 
ipo from {us)^. 



8.2. Momentum Transport 

To flnd an evolution equation for u = u}{tlj)R'^'V(j), we multiply (154) by nig {v ■ V0) R^, 
where v = w + u, integrate over w, and sum over species to flnd 



d 

— ^msUsR^ujiip) + 

s 

-^mslfd^w {5a, ■ V0) R^5f, 



V ■ (ris + msUsUsU + msUsuUt 

s 



I turb 



■(V0) R^ 
(jxS)-(V0)i?2 = ^(-\ 



;i73) 



where the viscous stress tensor is 

= Jd^WTHsWwFs, 

the source of toroidal angular momentum is 



n) 



(174) 



(175) 



where S'i"'' is given by (164) and we have used the quasineutrality condition (52). 
The second term in (175) is due to injected particles joining the mean flow, while 

% Explicit expressions for the collisional fluxes can be found in [58, 29] and in [36, 37, 35] for the case 
of non-zero u}{ip). 

+ Having determined the particle flux, we note that the classical, neoclassical and turbulent 
contributions are all independently ambipolar [I]. This is explicitly proved in Appendix C.2. 
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the first term contains also contributions from particle-conserving momentum injection 
meclianisms, e.g., momentum injection by waves [j!)]. 

We now wisli to write all terms except the time derivative and the source as full 
divergences. The second term on the left-hand side of (173) can be written so because, 
for any symmetric tensor field A, 

(V ■ A) ■ (V0) R^ = V ■ (i^^A ■ V0) . (176) 

The third term (due to the fluctuations) can also be written as a divergence: 



J^rrisQd^w (5a,- V0) R^6f, 



turb 



1 

C 

V- 



{6j X 6B),^^^ ■ (V0) R 



(177) 



6B' 



I 



(V0) R' 



turb 



Sir Att 

where I is the unit dyadic and we have used (156) for 6as, the fluctuating quasineutrality 
condition (146), expressed the Lorentz force as the divergence of the Maxwell stress in 
the usual way and used (176). The final term on the left-hand side can similarly be 
written as the divergence of the Maxwell stress associated with the mean magnetic field: 



--(jxB)-{V<P)R' 

c 



BB 



(V0) R' 



(178) 



_87r Att 

We now flux-surface average (173), using the identities (39) and (45) for the flux- 
surface average of a divergence and of a time derivative. This gives 



(179) 



where 



J = ^ms (usR^) 



(180) 



is the flux-surface-averaged moment of inertia and the total flux of angular momentum 
is given by 



IT. 



J2 W) ■ ■ (V0) R^ + J2 msOoWR'T, 



- - ■ (SBSB),^^^ . (V0) R' 



:i8ii 



where Tg is given by (167). 

Proceeding analogously to the derivation of the particle transport (see (168)), we 
can write the radial angular momentum flux due to the particles of species s as 



(W) ■ n, ■ (V0) R^ + msnsUj{i))R^Us ■Vip = j £wBms^ [R^v ■ Vcf)'' 
where v = u + w. We evaluate this explicitly in Appendix C.3 to flnd 



tot 



n 



^ \ ^EM I ^ 



(182) 
(183) 
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where ni^'^^ contains the classical, neoclassical and turbulent viscous stresses: 
= (^i^'^>^ + (^^% + (( [d'wmM'iiw . V0) KV^)^ . ) 

\ \ J I turb/ t/i 



with 

c 



and 



d^w'^w^w^CXF^,])^ : (V0) (V0) - ^ 



52 



d?wm 



J(^/')w|l f w X b 



B 



1^ iv^r 

m_, ;;7^ h /is 



52 



fi2 
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;i84) 



;i85) 



;i86) 



msUj{ip) {R^Ts)^ is the convective flux of angular momentum with {R?Ts)^ given by (cf. 
(170)): 



{R'T,)^ = (r' jd'w . VV^) C[Fo.] 



and the electromagnetic angular momentum flux tt^'^'' is 



TT 



(V0) 

EM 



(Vip) ■ ( -^SBSB + -SjSA) ■ (V0) 



;i87) 



(188) 



The two terms in the last expression are, respectively, the Maxwell stress and the 
advection of electromagnetic momentum by the particles* 



-Sj5A = J2 [d'w{hsw),^5A. 

C J c 



(189) 



8.3. Heat Transport and Heating 

Turning to heat transport next, we multiply (154) by rrisVLp' / 2 and integrate over all 
velocities to obtain 

3 (9 / 1 dA \ 

-— ra^Ts+V ■ Qs + ZsC ( Vv?o + ) ■ {"^sUs) + msUsU ■ (Vw) -Us + Us • Vu 



\j I turb 



(190) 



* In the long-wavelength limit, this second term is small as the fluctuating velocity is dominated by 
the E X B velocity, which carries no current. 
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where we have defined the heat fiux 

Qs = jd^w^rrisW^wFs, 
the colhsional energy transfer to (or from) species s 



and the energy source 



d^w-nisw'^Ss- 



(191) 



(192) 



(193) 



Acting in the same vein as for the particle and momentum transport equations, 
we now fiux-surface average (190). This flux-surface average is carried out in Appendix 
C.4, where we obtain: 

3 1 d 
2V' di 

Let us detail the terms in this equation. 



1 d 

V' (n,)^r, + — — V^' (g,)^ 



(194) 



8.3.1. The Heat Flux. The "heat flux" in (194) is deflned by the following formula, 
whose origin is explained in Appendix C.4. 6: 



turb- 



+ 



1 



-Ts + ZsCifo - -nisUj {i))R 



di)\ 



(195) 



This form of the heat flux is most useful for interpreting its constituent parts. The flrst 
term in the flrst line of (195) is the usual flux of thermal energy, the second term is the 
flux of the mean potential energy (electrostatic and centrifugal), and the third is the 
flux of the fluctuating electrostatic potential energy. The second line of (195) is minus 
the flow of heat and energy carried by the motion of the flux surfaces (remembering 
that ■ ViIj = —dip/dt and observing that the expression multiplying dip/dt is the 
enthalpy carried by the plasma^). As we include the potential energy fluxes in (195), qs 
is not, strictly speaking, the heat flux in the usual sense. We will elaborate upon this in 
Section 8.3.5 where we discuss the potential energy exchange term PP°*. For simplicity, 
we will continue to refer to {qs)^ as the radial heat flux. This usage will be vindicated 
on physical grounds by the appearence of {qs)^ multiplying the temperature gradient in 
the expression (235) for the entropy production. 

The first term is the enthalpy of the ideal gas of species s [ ] and the second and third terms are 
the potential energy contributions to the internal energy. 
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Unfortunately, (195) is not a useful form for actually evaluating the heat flux. An 
explicit form for it is calculated in terms of known quantities in Appendix C.5: 



+ 



jd?wes{h,V^)^-Vi) 



n^)^?^ (196) 

V> \-° lip 



turb / ip 

Let us now examine the heating terms on the right-hand side of (194). 



8.3.2. Viscous Heating. This is (see Appendix C.4.2) 

pvisc 



{r^^'% + mM^){R'rg)y^, (197) 

where (jii^'^''^ is given by (184) and (R'^Ts)^ by (187). Usually, the heating due to the 

mass flow (the second term inside the square brackets in (197)) is not included in the 
viscous heating. We group them together here because the term in the square brackets in 
(197) is precisely the species-dependent term in the total momentum flux (183). We will 
come back to this point when interpreting the rotational part of Pf°^ in Section 8.3.5. 
The combination of momentum fluxes seen in (197) will appear multiplying du/dip in 
the expression (235) for the entropy production, vindicating its interpretation as the 
total viscous heating of species s. 

8.3.3. Ohmic Heating. The mean Ohmic heating due to the induced parallel electric 
fleld is (see Appendix C.4.3) 

P^'^^ = K,{^){E.B)^, (198) 

where 

= ^ jd'ww\\F,,. (199) 

8.3.4. Compressional Heating. This is due to the motion of the flux surfaces and is 
given by (see Appendix C.4.4) 

Pr^" = -Ts (n.V ■ V^)^ = /n.V ■ (^^) ) , (200) 



dt |V^|^ 



where we have used the expression (44) for the velocity of a flux surface. 
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8. 3. 5. Exchange between Potential and Thermal Energy. The change in thermal energy 
due to exchange with the electrostatic and rotational potential energy is (see Appendix 
C.4.5) 

= - (^ZsC^o + ■ Vn,^ ^ - (Z,en,v?oV ■ V^)^ 

"^'^"^^ ^ V'ms{R'ns)^-lm,u\ilj){nsV^-VR')^ (201) 



+ 



2V' dt 



4> 

The first line is the energy exchange with the electrostatic potential energy, the second 
the exchange with the rotational energy. More transparently, we can group the terms 
in (201) involving (pQ together with the the terms involving ipo in our definition (195) of 
the heat fiux to write the total as 

^ ^ T/'/^ TT Y7 / ^ 7 (202) 

V'd^) \ dt I ^ 

= -Z,e {us (Us - V^) ■ Vv^o)^ , 

the work done by the mean electric field due to v^o-tt Unfortunately, we do not know the 
small (i.e., of order e'^Vths^s) poloidal component of UsUs; the only way of calculating 
the right-hand side of (202) is to calculate the left-hand side explicitly. This was the 
reason for splitting the work done into the "heating" and "fiux" terms and distributing 
them between (201) and (195), respectively. 

Similarly, gathering the terms involving the rotation in (201) and (195), as well as 
the second term in (197) (see discussion in Section 8.3.2), we can rewrite them as 

co\^) d_ 
2V' dt 



= ^m,ui'(i.){n,(U,-Vt)-VR^)^,. 

which is the work done by the centrifugal force (the Coriolis force does no work). Again, 
as with the work done by the electric field (202), we cannot explicitly evaluate the right- 
hand side of (203). 

ft Here, appears because we measure all velocities relative to V^p and so, for a force F, we consider 
UsF ■ to be work done accelerating the plasma to V^^, not heating. 
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8.3.6. Turbulent Heating. Finally, the energy exchange with the fluctuations is (see 
Appendix C.4.1) 



turb 



pdiss pdrive 

s s 



The last expression reflects the fact that P*"''*^ consists of two parts (see (E.28) in 
Appendix E.4): the turbulent heating (dissipation of fluctuations on collisions): 

Pf- = -(( ld?w^C[K]\ ) (205) 



\ \"' ^Os I turb/ ^ 

and the energy injection into the fluctuations due to the mean density, temperature and 
velocity gradients: 

odrive // fj3„../u \ V7„/,\ \ rr -^s 



J I turb / ^ 



dip 



+ (l^jd^wm.R^ (V0) ■ {vhsV^)^ ■ VV^^ 



turb / ip 

dui 

turb/ ^ 

This is the energy borrowed from the internal energy of the equilibrium and transferred 
into turbulence by such mechanisms as the ITG, PVG, ETG, and other gradient-driven 
instabilities [51, 61, 62, 63]. In Section 9, we will show that all of this energy is eventually 
returned to the equilibrium via the turbulent heating term: namely, after summation 
over species, net turbulent heating - the difference between pdnve P^^^^*^*^ - 

is either viscous heating (conversion of rotational energy into heat at constant internal 
energy) or due to turbulent eletromagnetic energy injected by antenna-like mechanisms, 
whose dissipation is the only piece of Ylis ^s^^^ ^'^^ cancelled by Ylis Pf™'^- 

8. 3. 7. Collisional Heating. Formally, the collisional heating term CI is both too large 
- 0{eVLsnsTs) or □(e^fi^'^s^s) if all mean temperatures are equal, as C[Pos] ~ e^^-^os if 
all temperatures are equal - and contains contributions from the gyrophase-independent 
part of F2S1 which we have not calculated. This problem has arisen from our decision 
to retain different temperatures for different species, a choice formally inconsistent with 
our orderings, as we acknowledged in the footnote preceding (85). There are two naive 
solitions to this problem. The flrst is to respect our ordering and assume that all 
temperatures are equal (which follows from (84)). In this case we can sum the heat 
transport equation (194) over all species and, as collisions conserve energy, the collisional 
heating term will vanish 

5^Cf) = 0. (207) 
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Alternatively, if, completely formally, we order the interspecies collision frequency to be 
small Usgi ~ e^z/ss, then this heating term will simply be the temperature equilibration 
between the Maxwellian equilibrium distributions of the species [43]: 

CP = E'^S(7^s'-7^.), (208) 

s' 

where 

(m^T^/ + rris'Ts) 

and InAg^/ is the Coulomb logarithm. However, this stretches the credibility of our 
ordering as there is no physical reason for i/gs' to be smaller than Ugg in general. Only 
in the special case of a very light species (electrons) or a very highly charged species 
(large- impurities) is it reasonable to so order the temperature equilibration time. 
Thus, we elect to treat electrons separately but sum (194) over all other thermal species 
(ions), which should all have the same temperature. 



9. Energy Conservation in Multiscale Gyrokinetics 

In this section, we bring together the results of Sections 7 and 8 to determine how energy 
and entropy flow between mean and fluctuating quantities. 

Firstly, in Section 9.1, we discuss energy conservation on the transport timescale. 
We show that the conserved quantity is the total energy of the plasma - the sum of the 
thermal energy of each species, the rotational energy and the magnetic energy. We also 
show that the fluctuations can do no net work upon the plasma, i.e., there is no net 
turbulent heating. | 

Secondly, in Section 9.2, we derive the conservation laws that hold on the fluctuation 
timescale. We show that the conserved quantity is the free energy of the fluctuations. 
We then show that its conservation is local, i.e., neighbouring flux surfaces do not 
exchange free energy. We also use the steady-state version of this conservation law to 
interpret the absence of work done by the fluctuations (demonstrated in Section 9.1) as a 
precise balance between two energy flows: kinetic energy of the thermal particle motion 
and of the bulk plasma is converted into turbulent fluctuation energy via instabilities 
driven by mean-field gradients (e.g., ITG [61, 62], ETG [51], PVG [63, 64, 65]); it is 
then cascaded to small scales in phase space, where it is dissipated (converted back into 
thermal energy) by collisions. 

Finally, in Section 9.3, we discuss why the gradients of the magnetic field cannot 
drive fluctuations on their own and note the impossibility of anomalous resistivity in 
our formalism. 

I This holds only in the absence of direct injection of energy into the fluctuating scales - as is (usually) 
the case in tokamaks. 
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9.1. Energy Conservation on the Transport Timescale 

9.1.1. Thermal Energy. We begin by considering the evolution of the total thermal 
energy. In Appendix E.l, we sum (194) over all species and simplify the result to 
obtain: 



3^ d_ 
2V' di 



V 




# 2V' dt 



+ (tt. 



du 



'EM 



/ ip dip 



J2 I ^"i^h 

s 

V'J+{E.j)^ 



(210) 



where J is the plasma's moment of inertia, defined in (180), Tr^o/'' is defined in (181), 
7T^^ in (188), and P*^''^ in (204). Examining the terms on the right-hand side of 
the above equation, we identify the first and second as exchange between thermal and 
rotational energy, the third as the Ohmic heating associated with the mean fields, and 
the last term as the energy injected by the heat and particle sources. We show in 
Appendix E.2 (and, by a different route, in Section 9.2.3) that the fourth and fifth terms 
on the right-hand side of (210), due to the fluctuations, cancel exactly (see (E.15)): 



turb 



+ 77. 



'EM 



/ ip dip 



0. 



(211) 



Thus, the heating due to turbulence arises entirely out of the turbulent contribution to 



TT. 



tot 5 



which results in viscous heating in (210) (the first term on the right-hand side) - 



converting rotational energy into thermal energy or vice versa. 



9.1.2. Rotational Energy. The last point is made manifest by considering the evolution 
of the rotational kinetic energy. Multiplying (179) by w(V'); "we obtain 



1 1 d 
2V' dt 



duj 



I i) dip 



V> dip 
a;'(V;) d_ 
2V' dt 



V'J + u{iP){S^-^)^ 



(212) 



where we have written the left-hand side as an exact time derivative and an exact spatial 
derivative and collected the terms arising from this manipulation on the right-hand side 
so as to parallel (210). We now add (212) to (210) to find 

1 d 

^ " \=^E.j), 

(213) 



V dt 



+ 



1 



5P +a;(^)5('^) - -m,u;2(V')i?'5i") 



Multiscale Gyrokinetics for Rotating Tokamak Plasmas 



41 



where the total kinetic energy of the plasma (thermal energy and rotational kinetic 
energy) is 

1 



U 



J2 9 ^""^^^ 



(214) 



and its flux is 




n 



dip\ 



^tot / 
/ ^ 



E 



Qs + -msUj'^{il))R^nsUs 



u 



(215) 



The two forms of (T^^^')_^ written above are shown to be equivalent in Appendix E.3. 
They serve two different purposes. The former is easier to evaluate for a specific 
practical implementation of the transport equations, while the latter is easier to interpret 
physically. In the second form of T^'^\ the terms are identified as: the heat flux (see 
footnote in Section 8.3 before (196)), the flux of rotational energy, the energy flux due 
to the viscous stress, the Poynting flux due to the fluctuations and a term due to the 
motion of flux surfaces. 



9.1.3. Magnetic Energy. The evolution equation (213) for U has two sources on the 
right-hand side: Ohmic heating and energy injection due to heat and particle sources. 
Let us now express the Ohmic heating term (the first term on the right-hand side of 
(213)) by using Poynting's theorem for the mean fields: 

Averaging this equation over a flux surface and adding it to (213), we obtain 



V' di 



msU^{'ijj)R'^Sl."'^ 

where we have used the flux-surface-average identities (39) and (45) and B x Vip = 
—I{iIj)B + R^B'^'Wcf) to rewrite the radial Poynting flux in the following way: 

(Ex B)-ViIj = -IN))E ■ B + R^B'^E ■ Vcj) = -I(ip)E ■ B + — (218) 

c ot 

Equation (217) is the desired energy conservation law on the transport timescale: there 
are no sources save those explicitly contained in the kinetic equation. 
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Note that, as (216) contains no contribution from the fluctuations, turbulence 
cannot cause any exchange between the total kinetic energy U and the energy of the 
mean magnetic field. These results, suprising at first glance, in fact follow directly from 
Poynting's theorem for the fluctuating fields: 

l^-f- + -^V ■ {6E X 6B) = -6E ■ 6j, (219) 
ot 8n An 

where we have dropped the energy in the electric field as it is negligible for non- 
relativistic plasmas. Averaging (219) over the fluctuations gives 

■ {6E X 6B\^^., = -{6E ■ 6j\^^^, (220) 

as the time derivative of the averaged fluctuating magnetic energy is smaller than the 
other two terms. This clearly demonstrates that the work done on the particles by the 
fluctuating flelds (the right-hand side) is balanced by the Poynting flux of fluctuating 
electromagnetic energy (the left-hand side), so the fluctuations cannot change the total 
kinetic energy U. The Poynting flux associated with the fluctuations appears explicitly 
in the second line of (215). If we wished to consider electromagnetic fluctuations driven 
by an external source, this would enter the energy- conservation equation not through 
a source term but through the fluctuating Poynting flux at the plasma boundary. We 
will return to the physical interpretation of this result in Section 9.3, after discussing 
the conservation of free energy in the next section. 



9.2. Free Energy Conservation and the Turbulent Cascade 

In this section, we derive the fluctuating counterpart to the energy conservation law of 
the previous section. The conserved (cascading) quantity for kinetic turbulence is not 
the energy of the fluctuations, but is in fact a combination of the entropy and the energy 
known as the free energy [11, 12, 13, 14]. In our notation, this quantity is 




(221) 



where we have averaged over an annular region of space centered on the flux surface 
labelled by ip (the composition of the flux-surface average and the perpendicular spatial 
average)!, but not over time. The first term in (221) is —Yls'^s^Sg, where 6Ss is the 
part of the mean entropy density of species s due to the fluctuations (note that this is 
not the same as the fluctuating part of the entropy density because 6Ss has a nonzero 
average; see Section 10.1). Thus, (221) agrees with the usual thermodynamic deflnition 
of Helmholtz free energy. 

I For the purposes of this section, we can visuaUse a flux surface as having a finite width (comparable 
to the intermediate length scale A introduced in Section 3.1) and so in our discussion we will drop the 
distinction between a flux surface and the annulus centered upon it. 
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In Appendix E.4, we derive an evolution equation for W from the gyrokinetic 
equation (144): 

d_ 
di 



d\nN, 3d\nTs 



dip 2 dip 



din To 



This is the free-energy balance on a given flux surface. 



—6B6B+-6A6j) 

4:71 C 



(222) 



du 
'(hp' 



9.2.1. The Local Cascade Law. Before interpreting the sources and sinks in this 
equation, we would like to highlight the fact that there is no flux of free energy - the 
turbulent cascade of free energy is essentially local and does not couple neighbouring 
flux surfaces. This means that turbulence excited on a given flux surface must dissipate 
on that surface and, moreover, that the turbulence is excited by the local gradients and 
cannot be excited by turbulence spreading from neighbouring flux surfaces. Succinctly: 
that which is stirred up in the flux surface, stays in the flux surface - the local (in space) 
cascade law of tokamak turbulence. 



9.2.2. Sources and Sinks of Free Energy. The first term on the right-hand side of (222) 
is the sink of free energy due to collisional dissipation§. The remaining source terms 
are the entropy generated by the turbulent fluxes transporting mean quantities along 
their gradients. As the fluxes are generally dominated by the largest turbulent scales 
and the collisional dissipation is most important when acting on the smallest scales (in 
phase space), we can interpret (222) in terms of a cascade in phase space analogous to 
the Kolmogorov cascade of energy in hydrodynamic turbulence [13, 1 1]: fluctuations are 
excited at large scales, then nonlinearly cascaded, whilst conserving W , to small scales, 
where they are destroyed by collisional dissipation. 

It is important to note that only gradients in the thermodynamic and mechanical 
quantities (density, temperature and angular velocity) appear in the right-hand side 
of (222). Only gradients in those quantities represent deviations from the global 
thermodynamic equilibrium and can, therefore, drive fluctuations (as the presence of 
fluctuations implies non-zero W). Gradients in the magnetic field strength or magnetic 
curvature do not have this property; the relaxation of the mean magnetic field to 
a minimum energy state proceeds independently of the fluctuations (i.e., there is no 
fluctuating contribution to (216); see further discussion in Section 9.3). 

§ The fact that it is a sink and not a source fohows from the i?- Theorem for the Unearised colhsion 
operator C[hs\ (see Section 10.2). 
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9.2.3. Statistically Steady-State Turbulence. Time averaging (222), we see that a 
balance between free-energy injection and dissipation is set up on a timescale much 
shorter than the intermediate time T (introduced for the purposes of time averaging in 
Section 3.1). This balance can be written as 

E^"'"-(4tf),^ = E^'"-. (223) 

s s 

where P^^^"™ and P/""" are defined by (206) and (205), respectively, and ir^^'^ by (188). 
This is another derivation of (211). This clearly demonstrates that the steady-state 
turbulence extracts all its energy from the equilibrium, but it also returns exactly 
the same amount of energy via collisions. This explains how the fluctuations can 
be continuously driven and dissipated without changing the total kinetic energy U. 
Turbulence can, however, move energy between various constituent parts of U, notably 
energy extracted from one species does not have to be returned to that species, and 
so fluctuations can redistribute energy between species (subject to the caveats about 
temperature equilibration at the end of Section 8.3). This can lead to systematic 
dependences of equilibrium temperatures on the nature of the species, e.g., in the case 
of effective heating of heavy minority ions [6G] . 

Similarly, turbulence can convert rotational energy into thermal energy (or even 
vice versa, e.g., if there is an angular momentum pinch or intrinsic rotation [67, 68]). 

Equation (222) implies that in order to sustain a non-zero steady-state flux of 
particles, heat or angular momentum, coUisional dissipation is required. This is the 
reason, alluded to in Section 3.4, why it was crucial to retain C[5/s] at the same order 
as d6fs/ dt. Similarly, this is why all simulations of the gyrokinetic equation (144) have 
to include some dissipation in order to reach a steady state [11, 69, 70]. 

9.3. The Special Status of the Mean Magnetic Field 

At several points in the above discussion, we have noted that the mean magnetic field 
behaves very differently from the other mean fields. In this section, we collect these 
individual points together to further elucidate the special role that the mean magnetic 
field plays in multiscale gyrokinetics. 

The mean magnetic field is determined from the Grad-Shafranov equation (136) 
and the evolution equation for the safety factor q^ip) (143) (see the detailed discussion 
in Section 7.3). These equations contain no direct contribution from the fluctuations. 
In this sense they are very different from the transport equations determining the mean 
density, velocity, and temperature profiles. This separation is also reflected in the 
energetics of the magnetic field. The fluctuations cannot directly convert magnetic 
energy into thermal energy (note the absence of any fluctuation-dependent term in 
(216)). Correspondingly, the magnetic field is not a source of turbulent free energy 
(there are no source terms in (222) arising from gradients in B), i.e., energy cannot be 
borrowed from the magnetic field to fuel fluctuations in the same way that it can be 
from the thermal or kinetic energy of the plasma. 
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Ultimately, these results imply that there is no turbulent contribution to the 
diffusion of mean magnetic flux. The fluctuations considered here cannot reconnect 
the mean magnetic field, i.e., they can change the topology of the total magnetic field 
B = B + 6B but only on small spatial and temporal scales, while the topology of 
B is preserved - indeed this is true automatically as we assumed axisymmetry of the 
equilibrium (see Section 3.2). Equivalently, we can say that there is no anomalous 
resistivity due to gyrokinetic turbulence (i.e., the evolution of q{ip), given by (143), 
cannot be sped up by turbulence). 

10. Entropy Flows and the Thermodynamics of Multiscale Gyrokinetics 

In the preceeding section, we have established the nature of the energy flows, both 
between the equilibrium fields and the fluctuations, and between neighbouring flux 
surfaces. In this section, we discuss these results in explicitly thermodynamic language. 

In Section 10.1, we derive an evolution equation for the mean entropy and show 
that this is not only consistent with the above results but also provides further 
evidence that flux surfaces interact only through particle, momentum and heat fluxes 
(as opposed to turbulence spreading). In Section 10.2, we relate the evolution of 
the mean entropy to notions of dissipation and equilibration through Boltzmann's H- 
Theorem. In Section 10.3, we relate the sources of mean entropy to sinks of free energy 
discussed previously. We discover that dissipation increases mean entropy via three 
routes: coUisional and turbulent temperature equilibration, Ohmic heating due to the 
induced electric field, and fluxes that strive to flatten gradients in the thermodynamic 
quantities n^, T, and u^ip)- We also discuss the effects that can occur in systems where 
two or more equilibration mechanisms compete. 

10.1. Entropy Balance 

The (Boltzmann) entropy density of species s is given by 



where h is Planck's constant. || To lowest order in e, fs = Fqs and so the mean entropy 
density of the plasma is 



where tt-q^ = {nisTs/ 2'Kh?) is the quantum density of states for the ideal gas 
of species s. This is just the standard result for a mixture of ideal gases in 
local thermodynamic equilibrium [(iO], as expected. We are in local thermodynamic 
equilibrium: fg is Maxwellian to lowest order (the equilibrium is local, not global, as 
rig and T, are functions of space). Because we are considering a given flux surface, the 

II For an explanation of this suprise appearance of /i in a classical formula see §7 of [60]. 




s 



(224) 




(225) 
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presence of the rigid rotation u{tp) does not contribute to the entropy density at this 
surface. Expanding fs = Fg + 6fs, we find an expression for the mean entropy density 
that retains terms up to order e^: 



H 



E 



d wF^ In 



STT^h^F, 



mi 



2F, 



(226) 



turb 



where the first term is the entropy of the mean distribution Fs and the second term is 
the contribution from the fluctuations, 5Ss, which was introduced in the discussion of 
free-energy conservation (see (221); in (226), 5Ss is also time-averaged). 

We now derive an evolution equation for H. Multiplying the kinetic equation (1) by 
— [1 + In [Sn^h^fs/ fn^)], integrating over all velocities, summing over species, averaging 
over the fluctuations and over a flux surface, we obtain 
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V' dt 
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where we have used (45) to simplify the time derivative, (39) for the flux-surface average 
of a divergence, and defined the radial entropy flux 
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the coUisional entropy production 
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(230) 



In Appendix E.5, we show that the entropy flux (228) can be written in terms of 
the particle and heat fluxes through a flux surface as follows 



Jr'a^-E 



(231) 



where we have introduced the chemical potential of an ideal gas of particles of species s: 



T,(^)=T,ln 



Qs 



TJn 



«Q6 



(232) 



(we used (96) to obtain the last equality; note that the last two terms of the first 
equality are the potential energy per particle). Equation (231) is analogous to the same 
result for a mixture of ideal gases [71, 72]. We have already seen, in Section 8, that each 
annulus interacts with neighbouring annuli via the particle, momentum, and heat fluxes. 
The expression (231) for the entropy flux shows that the annuli exchange entropy only 
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through the particle and heat fluxes (but not momentum, as we should expect, given 
that the mean entropy density does not depend upon u{tp)). 

By expanding fg in terms of Fg and Sfs, we can write the mean (coUisional) entropy 
production term as 



where we have used (104) for Sfs- In this expression, we can clearly identify the separate 
contributions from collisional dissipation of Fg and coUisional dissipation of 6fs. It 
is important to note that the second term in (233) is minus the time average of the 
first term on the right hand side of (222) - the rate of entropy production due to the 
fluctuations is precisely the rate at which free energy is dissipated. 

We determined in Section 9 that the energy budget for the fluctuations consists of 
borrowing some energy from the mean thermal or bulk kinetic energy, cascading it to 
small scales, and then returning it. We also showed that in this process no net work was 
done on the plasma. We seem to have arrived at a contradictory result: we know that, 
in the absence of viscous heating, there is no local heating, but (233) shows that we have 
a monotonically increasing entropy (the right-hand side of (233) is positive-definite if 
we are away from thermodynamic equilibrium). This apparent contradiction is resolved 
in the next two sections, where we consider the precise nature of the thermodynamic 
equilibration of our system. 

10.2. Boltzmann's H- Theorem and Global Equilibrium 

Boltzmann's if-Theorem [73] states that for a homogenous dilute gas of hard spheres, the 
entropy H^, as defined by (224), is non-decreasing and the only steady-state distribution 
is a Maxwellian. Let us examine the corresponding theorem for the plasmas under 
consideration here. The general form of the if-Theorem for a dilute plasma is: if 
fs{r,v) is the distribution function for species s, then [43] 



with equality achieved if and only if for each species fs is a Maxwellian distribution and 
all their temperatures and mean velocities are equal. We have already called upon this 
result in Section 6.1 to prove that Fqs is Maxwellian. We now examine its implications 
for the thermodynamics of the plasma. 

First of all, we see that C^^^ > 0, i.e., collisions generate entropy. Moreover, the two 
terms in the expression (233) for C^^^ are separately non-negative. Thus, the relaxation 
of the equilibrium Fg and the dissipation of the turbulent free energy are separate net 
sources of entropy. 




(233) 




(234) 
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It is instructive to discuss the equilibrium towards which these relaxation process 
are driving the system. In the absence of sources or fluxes of entropy connecting the 
plasma to the outside world, the system is in a global equilibrium if C*^^-* vanishes when 
integrated over the plasma volume. However, as {C^^^)_^ is positive definite on each flux 
surface, for its volume integral to vanish, it must vanish locally on each surface. From 
the if-Theorem, we know that this is the case only if fg is a Maxwellian up to possible 
errors of 0{e^fs) (these smaller deviations from a Maxwellian generate an evolution 
of the entropy on timescales longer than te and so are negligible). By applying the 
i^-theorem to the first term of (233), we see that Fg must be Maxwellian for C^^^ to 
vanish. From the explicit expression (B.7) for (Appendix B), we conclude that for 
to be Maxwellian, all the gradients in F^g (i.e. those of Ug, Tg and u) must vanish and 
the neoclassical distribution function f]"^^^ must also be Maxwellian (in the absence of 
gradients in Fog, (119), which determines Fs^^\ admits purely Maxwellian solutions).^ 
Thus, collisions, which drove the distribution on each flux surface to be Maxwellian on 
shorter timescales, now strive to smooth out the gradients between the flux surfaces. 

The if-theorem has determined the form of the background distribution function, 
the direction of the free-energy cascade from injection scales to coUisional scales in phase 
space, and it has determined the ultimate equilibrium to which the system wishes, in 
some sense, to arrive at. There remains but one issue to clarify: how the relaxation 
towards this global equilibrium is achieved. 



10.3. Approach to Equilibrium and Multichannel Transport. 



The precise nature of the relaxation process is hidden inside the entropy production 
terms in (233). In Appendix E.6, we show that the collisional dissipation in (233) 
can be expressed entirely in terms of fluxes through a flux surface and local relaxation 
termsl [1, 2]: 
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P,"*^- in (198), (F,)^ in (170), T, in (232), 
in (184), and (P^F,)^ in (187). Note that 
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We now interpret (235), term by term. 



(236) 



% Even though we do not have an evolution equation for it, f]"'^'' evolves in time due to the time 
evolution of the right-hand side of (119) and the time evolution of the collision frequency. 
I We note that this is a somewhat arbitrary distinction. For example, the viscous heating (the last 
term in (235)) can be considered as a flux-gradient term for the flux of angular momentum, or a local 
heating term. 
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The first term is just the collisional equilibration of mean temperatures; it vanishes 
if all temperatures are equal. 

The second term contains two effects. Firstly, if all temperatures are equal 
{Ts = T), then we can use (211) to write this term as a viscous heating term due 
to ( 

^EM / combine it with the last term in (235) to form the total viscous heating 

— (1/T) (jittf''^ du/ dip. Secondly, if there is no rotation (w = 0) then we find that, 

from (211), -Pg"'^'^ = 0. Therefore, in the absence of rotation, this term represents a 
change in entropy not due to net heating but due to the differential heating of different 
species by the turbulence. Hence, in the absence of rotation, this term is a source of 
entropy if it equilibrates temperatures but a sink of entropy if it drives them apart. 
That this term can generate entropy without net heating is part of the resolution of the 
apparent contradiction noted at the end of Section 10.1. 

The third term in (235) is the Ohmic heating due to the induced parallel electric 

field. 

The remaining terms in (235) (the second sum) can be interpreted as fluxes 
multiplying their corresponding gradients. 

The first term contains the particle flux and the gradient of the local chemical 
potential. The appearance of the chemical potential is natural as we can consider the 
flux surface to be a thermodynamic system that can gain or lose particles. 

The second term is just the heat flux (in the extended definition adopted in 
Section 8.3.1) multiplying the temperature gradient. 

Finally, the term multiplying the angular velocity gradient is not precisely the 
momentum flux as the electromagnetic momentum flux vr^^^ is hidden within P*"'''^ as 
discussed above. This ambiguity arises because, if the temperatures are unequal, it is 
impossible to separate uniquely the temperature equilibration due to turbulence driven 
by angular velocity shear and viscous heating due to the angular momentum flux driven 
by the turbulence. 

Individually, these flux-gradient terms are positive for any flux that acts to flatten 
the corresponding gradient, and negative if the flux steepens the gradient. | The presence 
of turbulent contributions to the fluxes in these terms is the second part of the resolution 
of the apparent contradiction that was raised at the end of Section 10.1 - the fluctuations 
can increase entropy without heating the plasma by flattening gradients. From (234), 
we know that C^^^ is positive. However, it is clearly possible to have some negative 
terms in (235) without violating this constraint - up-gradient fluxes of particles, heat, 
or momentum are examples of just such negative terms. It is the multi-channel nature 

X Take, for example, the first term in the sum, which is the entropy production arising from the 
particle flux (Fg)^ and the density gradient din Ng/ dip; assume all other gradients are zero. If tp 
increases away from the magnetic axis, (r^)^ will be positive for a flux or particles that is going away 
from the magnetic axis and dlnNg/ dip will be negative for a density proflle that is peaked on-axis. 
Thus, their contribution to the entropy production — X^s (^s)^ (din Ng/ dip) will be positive for a flux 
that relaxes the gradient, and negative for one that steepens it. This is independent of the definition 
of ip as it enters in both terms and the effect of flipping the sign of ip cancels between the two. 
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of the transport that allows such phenomena: offsetting a pinch in one channel against 
a larger outward flux in another to observe a net increase in entropy. 

Stating that the transport is mult i- channel simply means that gradients in any 
one thermodynamic quantity (n^, Tg, or u) can drive fluxes of any other quantity. 
For example, trapped-electron-mode turbulence can be driven predominantly by the 
electron density gradient but produces both a particle flux and an electron heat 
flux [ ]. Combining this multi-channel property and the fact that the terms due to 
the fluctuations in (235) can have differing signs and are only constrained to be positive 
upon summation can result in suprising effects - such as the spontaneous steepening of 
a gradient in u merely from temperature-gradient driven turbulence [Gy]. 

Finally, we note that, as expected from the discussions in Section 9.3, gradients of 
the mean magnetic field do not appear as sources of entropy. 

11. Multiscale Gyrokinetics at Low Mach Number 

In Sections 4-8, a multiscale hierarchy of equations is derived from the Fokker-Planck 
kinetic equation by a systematic expansion in e = p/a. The full generality of these 
equations is, in fact, unnecessary for a large class of plasma conditions, namely those 
where the Mach number M = Ruj{ip)/ Vths of the toroidal rotation is low, i.e., M ^ 1. 
In this section, we investigate the low-Mach- number limit: formulating the expansion 
in M ^ 1 precisely in Section 11.1 before presenting the results of this expansion in 
Sections 11.2-11.6. It is not necessary to read the detailed derivations in Sections 4- 
9 and their attendant appendices to use the results presented this section; however 
familiarity with the notation of Sections 2 and 3 will be assumed. Because of the 
considerable simplifications that the low-Mach-number expansion brings, the exposition 
in this section can be read as an easier-to-grasp summary of the basic structure of 
multiscale gyrokinetics. 

11.1. The Low-Mach-Number Ordering 

There are two distinct low-Mach-number regimes. In the extreme limit of M ~ e, all 
dependence on the angular velocity u drops out of the equations for both the mean 
and fluctuating fields, and (179) is no longer sufficient to calculate angular momentum 
transport. This is the "low-flow regime" discussed in detail in [75, 25]. In this paper, we 
deal instead with the intermediate limit where e <^ M <^ 1. This is the expected 
parameter regime for many current and future fusion experiments (see Table 1 in 
Section 1). 

There are two physical effects of plasma rotation that survive in this limit. They are, 
firstly, the suppression of small-scale turbulence by the perpendicular flow shear [76, 77] $ 

I The physics of the suppression of turbulence by sheared flows was originaUy discussed in these papers 
in the context of poloidal flows. The crucial fact is not that the flow shear is poloidal, but that it is 
perpendicular to the magnetic field. Thus, the same physical mechanisms (shearing apart of eddies, 
tilting of eddies due to the flow shear) apply here, despite the fact that the perpendicular shear comes 
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(the u ■ dhs/ dRs term in (144)) and, secondly, the instabihty drive due to the parallel 
flow shear that can itself excite turbulence [63, 64, 65] (the du/dip term on the right- 
hand-side of (144)). We retain the first of these terms by ordering the fluctuating 
frequency in the plasma frame 

|.„.V)^^ ,237, 

independently of the Mach number. The shear in the flow is hidden within the spatial 
dependence of the velocity. To make the presence of perpendicular flow shear more 
transparent, we assume that the domain of interest is shorter than the length scale of 
the flow and then expand the velocity around its value uq at the centre of the domain 

RsO- 

d d \ dh ^'''^ 

Ft^'''' OR J ^ ~ ■ [^""^^^"^l ■ + ■ ■ ■ ■ 

The second term on the right-hand side of this expression is now explicity the action 
of the perpendicular flow shear upon the distribution function hg. To retain only the 
part of this term due to shear in cj, we formally order the gradients of w ('?/') to be sharp: 
aV In u:{il)) ~ so that du/ dip is zeroth-order in M. For conciseness, we will continue 
to write the advection term unexpanded as u{Rs) d/ dRs, but the reader should keep 
in mind the idea of solving these equations in a small domain and Taylor-expanding the 
mean velocity field. 

This ordering for the gradients of u is also precisely the ordering required to 
retain the parallel-velocity-gradient drive term on the right-hand side of the gyrokinetic 
equation. We also order the timescale of variation of uj{ip,t) as Mte rather than r^, in 
order to retain the transport-time variation of the velocity. With these orderings, we 
will be able to keep the two effects we wish to examine and neglect all other effects of 
rapid toroidal rotation, e.g., Coriolis and centrifugal forces. In the following sections, 
we will apply this low-Mach-number expansion to the multiscale hierarchy derived in 
the preceding sections and present the ensuing equations to leading order in M. 

11.2. Mean Fields 

First we derive the equations for the mean distibution function Fg and the mean fields 
E and B. 

In the low-Mach-number limit, the mean distribution function is given by (see 
Section 11.3 for details) 

Fs = Fos iHRs),es) + Fis (i?„£„/i„a) + + O (eM/,) , (239) 



Fos = ns{ij{Rs)) 



rris 



_27rT,(^(i?,)) 
from a purely toroidal flow. 



3/2 



^-es/nWRs))^ (240) 
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where Ug and Tg are the density and temperature of species s, the energy variable is 
now given by (cf. (70)) 

= ^mV + 0(M2T,), (241) 

and Rg and fig are still given by (64) and (66), respectively. The distinction between 
the particle velocity v and the peculiar velocity w = v — u can now be dropped, and so 
we use w as the particle velocity throughout this section. The neoclassical part of the 
mean distribution function. Fig, is found from (118) and (120)-(122), as before. 
The mean electric field is 

E = -V$(^) + O (eM^^B) , (242) 



c 

where $ is related to the angular velocity (^{ip) through (61). The mean magnetic field 
is given by the two functions ip{R, z) and via (31). In the low-Mach-number limit, 
(136) becomes the usual Grad-Shafranov equation: 

(243) 

where pg = UgTg is the pressure of species s, which is now a flux function (see 
Section 11.3). The second component of the magnetic field, /(V^), is determined by 
the same method as in the high-flow regime - this is detailed in Section 7.3. 



11.3. Vanishing of Centrifugal Effects 

As an example of how these results have been obtained, let us prove explicitly that ipo 
and the poloidal variation of the density (see Section 6.3) vanish in the low-Mach-number 
limit. First, we expand the lowest-order quasineutrality condition (97) in powers of M 
to find 

'^ZgNg{ip)exp 

s 

Since the only spatial dependence in (244) is via tp, we know that the solution (for (po, 
given Ng) will be of the form = + 0{TgM'^ / Zgc). However, we have defined 

V?o so that it has no fiux-surface average (see (58) and (59)). Thus, ipQ must vanish to 
lowest order in M^. Therefore, from (96), we find rig = A^s(V') + O(M^ns). Thus, Ug is 
a fiux function to the required order in M and we can drop the distinction between A^^ 
and Ug. The quasineutrality constraint (244) then reduces to 

^Z,n,(^)= 0. (245) 

s 

11.4. Fluctuations 

In the low-Mach-number limit, the fiuctuating distribution function is (cf. (104)) 

6fg = -^6^{r)Fog + hg {Rg, Eg, fig, a) + O (eM/) , (246) 



Zg&fo 
Tgi^) 



O {M^Ug) 



(244) 
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where 6ip is the fluctuating electrostatic potential, which determines the fluctuating 
electric field to the required order: 



5E = -V5ip + O i^eMB 



(247) 



The gyrokinetic distribution function hs is given by the low-Mach-number gyrokinetic 
equation, which is the same as (144) but with centrifugal and Coriolis effects neglected: 



d_ 
di 



d 



OR. 



dhs 
OR.. 



{C[h,]) 



R 



d_ 
di 



+ u(Rs 



d 



OR. 



R 



dFos ^ mJ{ilj)w\iFos dco 



dip 



dip 



(248) 



where the guiding-centre drift velocity is now (cf. (Ill)) 

= 7^ X {wlb-Vh+-w\V\YiB\+0{eMvth^ 



The definition of x is now 

1 



-w ■ 6A, 



c 



(249) 



(250) 



because Sip' ~ 6(f in the low-Mach-number limit. The definitions (110) of and (72) 
of the gyroaverage remain unchanged. 

The equation for hg is closed through constitutive relations for the fluctuating fields. 
The fluctutating potential 6(p obeys the quasineutrality condition (cf. (146)): 



E 



(251) 



and the fluctuating magnetic field is 6B = 5B\\h + b x V(5v4||, with 5A\\ and 5B\\ 
determined from the parallel (149) and perpendicular (153) components of Ampere's 
law (153), respectively. 



11.5. Transport 



The long-time evolution of the densities Ug and temperatures Tg are given by the low- 
Mach-number limits of the transport equations (166), (179), and (194) of Section 8:§ 

l_ d_ 
V' di 



1. ^ 
V' di 

3^ d_ 
W' dt 



c(n)\ 



(252) 
(253) 



pvisc ^ pOhm ^ pcomp ^ pturb ^ (Cf^^) + (^i''^ • 



(254) 



§ With the transport-time variation of ijj{4>) ordered as M ^ as discussed in Section 11.1, aU terms in 
(179) are the same order. 
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The expressions (170) for the particle flux (F^)^, (164) for the particle source (180) 
for the moment of inertia J, and (193) for the energy source Ss are all unchanged in 
the low-Mach- number limit. The momentum flux (183), source (175), and the heat flux 
(196) lose some terms that are higher-order in M: 



<T}=Ei-i'% + {<'),^ (255) 



where (^iri^'^'^j and n^j^^ are given by (184) and (188), respectively, 
5^^) =J2[ d^'^^s {w ■ V0) R^Ss, (256) 

jd'wss {^^^ ■ Vi^ C[Fo,]y + (^Jd^wSsF^^'Wns ■ VV^^ (257) 



5 n^T(l^^^^^^^ 



J I turb/ i/) 



^ \^ I il, \ \ J I turb/ V 

The Ohmic heating pO"^™ is given by (198), the compressional heating P™™p by (200); 
they are unchanged in the low-Mach-number limit. The other two heating terms, viz., 
viscous (197) and turbulent (204), are simplifled: 

duj 
dip' 



Pr = -{-i^%^, (258) 



-Z,e(( ld^w(h, + ■ V ) X 

pAiss ^drive 



r/ turb/ ^ (259) 



where (^ni^^^^ is given by (184) and P^'^^^ and P^'^"™ are given by (205) and (206), 
respectively. 

The system of equations (252)-(259) with the solutions for the neoclassical 
(Section 11.2) and fluctuating (Section 11.4) distribution functions and the mean 
magnetic fleld (Section 11.2) is a closed system for coupled turbulence and transport in 
the low-Mach-number limit. 



11.6. Energy and Entropy 

In the low-Mach-number limit, we are able to neglect the rotational energy compared 
to the thermal energy of the system (despite the viscous heating appearing in the 
temperature evolution equation (254) - because we allow the turbulence to generate a 
large momentum flux from the sharp gradients in u{'iIj)). Thus, the energy conservation 
law (217) reads 



V' di 




(260) 
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where the flux of kinetic energy (in the low-Mach- number limit, this is purely thermal 
energy), previously given by (215), is 




(261) 



Here Qg is the usual heat flux deflned by (191) and (qs)^ is given by (257). 

The free energy (221) of the fluctuations continues to satisfy (222) in the low- 
Mach-number limit, but with all instances of Ng replaced by Ug. Note that because the 
gradients of uj{tp) have been retained in the low-Mach-number ordering, they contribute 
to the free-energy balance, and hence can drive turbulence. Equations (260) and (222) 
should be considered as constraints on the low-Mach-number system of Sections 11.2- 
11.5: any solution of these equations should identically satisfy these constraints. 

The entropy balance of Section 10 is unaffected by the low-Mach-number expansion 
(up to changes in the chemical potential (232) and the transport fluxes). Importantly, 
the entropy generation and source of free energy due to the angular velocity gradient 
remains. This is understandable as we have constructed our ordering so that the 
momentum flux is flnite and the gradient in the flow is flnite, and thus the entropy 
generation due to viscous heating remains flnite and comparable to all other forms of 
entropy generation. 



12. Summary 

In this paper, we have presented a complete and self-consistent derivation of multiscale 
gyrokinetics from flrst principles, drawing on and in some cases correcting previous 
work [1, 2, 3, 4, 5, 6, 38], as well as incorporating the evolution of the equilibrium 
magnetic fleld into this theoretical framework. We have shown how, in this framework, 
the gyrokinetic equation for fluctuations can be derived [2, 10, 38, 78] (Section 7.4), in 
parallel with the neoclassical drift kinetic equation (Section 7.1) more usually derived 
in a different setting [35]. This system of equations describes fast, small-scale and slow, 
large-scale perturbations of the local Maxwellian equilibrium driven by the gradients 
in the local equilibrium. We close this system on the long (transport) time scale 
by derving evolution equations for the equilibrium magnetic field (Section 7.3) and 
equations governing the transport of particles (Section 8.1), momentum (Section 8.2) 
and energy (Section 8.3). 

We are now in a position to draw together all of these results to form a clear picture 
of the evolution of a rotating turbulent plasma in a tokamak. The summary of our work 
can be reduced to three key points. 

Firstly, flux surfaces form effectively isolated systems on timescales shorter than 
the transport time. They rotate toroidally as rigid bodies (62), they are in local 
thermodynamic equilibrium (93) and the turbulence within them is established via a 
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local balance between energy injection and dissipation (223). Thus, we can conceptually 
think of a tokamak plasma as being made up of annular regions with a radial width given 
by some intermediate length (14) - these regions are independent on the intermediate 
time scale (16). It is this decomposition that underlies analytical and numerical 
considerations of the gyrokinetic system (144), (146), (149), and (153) in any one such 
annular region, the so-called "local approximation" [79]. 

Secondly, on long timescales (48), these independent annular regions are coupled by 
the transport equations of Section 8. These are local equations for the plasma density 
(166), angular velocity (179), and temperature (194), linking each annular region to 
its neighbours via transport fluxes that can be computed from averages over local 
solutions to the small-scale equations. Despite the dramatic simplification that has 
occured in going from the full kinetic equation (1) to evolving only a few flux functions, 
we have retained all the critical information about the plasma. This linked-fiux-surface 
approach to transport provides both a transparent theoretical framework and a basis 
for practially feasible numerical simulations of global tokamak dynamics on confinement 
timescales [6, 23, 24]. 

Thirdly, the (global) mean fields and the (local) fluctuations are also linked via 
energy and entropy flows (and so are individual flux surfaces). Section 9 is dedicated to 
the flows of energy, which express the basic conservation properties of the system. The 
results of this section are consequences of the multiscale representation of the plasma 
dynamics that is derived in the preceding sections. No extra information was used 
in their derivation. Thus, the energy balance (217) is a statement that if the mean 
quantities are evolved according to the transport equations of Section 8, then the 
mean energy of the system is conserved. Similarly, if the gyrokinetic equation (144) 
and the field equations (146), (149), and (153) are solved to determine the fluctuating 
distribution function and fluctuating electromagnetic fields, then the free energy (221) 
is the conserved quantity in the sense that (222) is satisfied by the solutions. These two 
statements are what is meant by energy conservation in our multiscale system. 

A direct consequence of the multiscale plasma dynamics is the evolution of our 
multiscale system towards thermodynamic equilibrium - the entropy of the system 
evolves according to (227). As a consequence of writing the entropy production in a 
convenient form (235), we are able to diagnose how our system proceeds towards a global 
thermodynamic equilibrium and what that equilibrium is. As the system approaches 
a global equilibrium, the coUisional and (relative) turbulent heating equilibrate the 
temperatures and the Ohmic heating dissipates the parallel induced electric field; the 
fluxes flatten gradients. Therefore, the global equilibrium is a state with no temperature 
differences, no parallel induced electric field, and no spatial gradients in n^, Tg or u 
(as a consequence of this, there are no fluctuations). However, this simple diagnosis 
is complicated by the multichannel nature of the transport in our system: far from 
equilibrium, it may be thermodynamically favourable for the system to initially steepen 
gradients in some of Ug, Tg or u in order to flatten other gradients more rapidly. 

In the course of our presentation of this general framework, we have elucidated a 
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number of issues that have the potential to be confusing and indeed have been so in the 
past. It is useful to itemize some of them here: 

• Turbulence is local to a given flux surface and does not spread - Section 9.2. 

• There is no turbulent heating - Section 9.1.1 (except in the sense of turbulent 
redistribution of energy between species). 

• The evolution of the mean magnetic field is determined only by the mean profiles, 
and not directly by the turbulence - Section 7.3. 

• Gradients in the magnetic field do not drive fluctuations - Section 9.3 

• The absence of the parallel nonlinearity in the gyrokinetic equation (144) does not 
upset the energy or entropy conservation properties of the system - see footnote on 
page 27. 

Let us now discuss some further directions and natural extensions following from 
this work. 

We have accomplished the construction of multiscale gyrokinetics by exploiting 
the scale separations inherent to plasma turbulence: between the fluctuations and the 
equilibrium and between the particle motion and the fluctuations (Section 3). However, 
these are not the only scale separations available to us. In [18], this framework is 
extended by exploiting the scale separation between ion- and electron-scale fluctuations, 
where this scale separation is found to have profound implications for the structure of 
the fluctuating magnetic field. 

In this paper, we have restricted ourselves to situations where the equilibrium 
distribution functions for all species are Maxwellian. In a burning plasma, we expect 
there to be at least one non-Maxwellian species - the fusion-produced a-particles. In 
[33], we consider precisely this situation - demonstrating how such a species can be 
incorporated into the theoretical framework developed here. 

We have also restricted ourselves to the so-called "high-flow ordering," where the 
mean rotation velocity is ordered to be comparable to the sound speed. The "low-flow" 
ordering, where the mean velocity is comparable to the drift velocities, requires a much 
more involved calculation to determine the transport of toroidal angular momentum 
correctly [25, 75, 80]. It is precisely this regime that must be considered if a detailed 
and complete study of the so-called "intrinsic" rotation is to be undertaken [81]. This 
alternative ordering may hold the key to how a slowly-moving or stationary plasma 
can transition into a rotating one [n2] or how L-H transitions occur in the edge of a 
slowly-rotating plasma [n>>]. 

Finally, it is impossible to present a set of equations such as the ones found in 
this paper without discussing how to solve them numerically. Currently, the linked- 
flux-tube code TRINITY [23] solves the low-Mach-number system of equations presented 
in Section 11. It is anticipated that this code will be extended to solve the complete 
set of equations derived here. With the advent of ever faster supercomputers, the 
combination of multiscale theoretical approaches and advanced numerical algorithms 
may allow us, for the first time, to simulate a burning plasma, from the small scales 
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of kinetic turbulence, through the scales associated with energetic fusion products, to 
the plasma equilibrium scales and even the timescales of the resistive evolution of the 
plasma. As a tool to study potential reactor designs, this would be invaluable. Hitherto, 
experiments have been optimised for MHD stability. Henceforth, it is possible to envision 
a time when they will be optimised not only for macro-stability but also for micro- 
stability and designed to achieve dramatically reduced turbulent transport (e.g., [84]). 
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Appendix A. Derivation of Rs, fis-, £s and i) 

Here we calculate the time derivatives of the variables Rg , /ig , and along the particle 
trajectory, which allows the derivation of the second order gyrokinetic equation (108) 
from the general kinetic equation (71). Note that we only require these derivatives 
correct to first order in e, because the higher-order (transport-timescale) effects are 
treated separately in Section 8. 

Appendix A.l. Some useful properties of the gyroaverage 

In performing gyroaverages in this Appendix, we will need to replace temporal and 
spatial derivatives at constant v with derivatives with respect to Rs, Eg, fis and i). For 
any quantity g{r, v) that has small-scale spatial structure (i.e. spatial variation on scales 



Ps), 



(yRs) ■ 
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+ (Ve.) 



+ (V/i,) 
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(A.l) 
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g + Oiek^g). 
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Similarly, for time derivatives, 



di 



d 



9 



dt 



g + Oie'^sg). (A.2) 

We will also use the fact that the gyroaverage commutes with derivatives in the 
Eg and Us variables. 

Another general result we will use is that 



dg_ 



Or 



Vg + 



Rs ,Ss ,f^s 



dv 

dg 



dg^ 
dv 



w±-Vg + (w X b) ■ h O(e^lsg) 

ow 



(A.3) 



w±_ ■ Vg + Vl. 



dg_ 
d'd 



OieQsg), 



where we have used the defintions (64), (65), (66) and (68) of Rg, Eg, Hs and to, 
respectively. Gyroaveraging the identity (A.3) at constant Rs or at constant r, we 
find that 



dg_ 



{w± ■ Vg)^ = Qs 



dg_ 



(A.4) 



R3,£s,^''3 



Finally, we will need the following identities in order to evaluate gyroaverages 
explicitly: 

{w^)^ = 0(et;thJ, {w^w^)^ =^(^\^bb) + 0(e<J, (A.5) 

where we have used (68) to carry out the integration over Similarly, we can use the 
definition (64) of Rg in terms of r to find, for any mean quantity g{r,Es, fJ^s), 



{g{r,Es,fis))R= {g{Rs,es,fis))R + (^-Yj— -^aj + 0(e^ 
= g{Rs,Es,fis) + 0{e^). 



(A.6) 



Appendix A.2. Derivation of Rg 

Rs is defined by (64). Taking the full time derivative along the particle trajectory, we 
get 

b X w 

I V 



Rs 



d ( b X w 

— r 

dt \ 



b X w . 
r V -V 



(A.7) 



where we have neglected terms of order 0(e'^fthj (time derivatives of the mean quantities 
Qs and b) and the spatial derivatives are taken at constant v. We now use r = v and 

w = i) — V ■ Vu = as + Stts — V ■ Vu, (A. 8) 
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where the accelerations are given by (155) and (156), to find 

c 1 
Rs = wub + u + —bxV{ipo + x) + -Bbx {w^ ■ V6A) 

+ — \h X {u ■ Vw) + b X (w ■ Vw) — {u ■ Vb) X It; — {w ■ Vb) x w 

(A.9) 

+ {bxw)w-V\YiB] + 0(e2t;thJ, 

where we have used u ■ Vlni? = (axisymmetry). The lowest-order component of this 
expression is 

Rs = w\\b + u + O {evtij . (A. 10) 

We now proceed to gyroaverage (A.9) term by term. We use (A. 6) on the first two 
terms to find| 

{w\\h + u)^ = w\\h{Rs) + u{R,) + 0(eVhJ. (A.ll) 

This is the combination of parallel streaming along the field line and convection by the 
mean fiow (but it is now the guiding centres, not the particles, that stream along the 
field line and are advected by the fiow). The third term in (A.9) gyroaverages to 

-|6 X V (<^o + X))^ = |b X V (<^o + (x) J = 1^ X V<^o + {V^)r. (A.12) 

where is defined by (110). These terms are the E x B drifts due to ipo and due to 
the averaged fiuctuating potential {x) r seen by the particle (x is just the electrostatic 
potential in the particle frame). The fourth term of (A.9) vanishes upon gyroaveraging: 

\-bx {w^-WSA)) =^bx {w^-W5A)j^ = Q, (A.13) 

B I R ^ 

where we have used (A. 4). The fifth term in (A.9) becomes the centrifugal drift: 

^hx{u- Vn)\ = -^h X {u?{^)RVR\ + 0(e2t;thJ, (A.14) 
where we have used u = uj{il))R^'W ip and expressed Vit as 

Wu = ^R^ (VV^) (V0) + u{^)R [(Vi?) (V0) - (V0) (VR)] . (A.15) 

(Jj iy 

Combining the gyroaverages of the sixth and seventh terms of (A.9), we find that they 
give rise to the Coriolis drift: 

[6 X (it; ■ Vu) — {u ■ Vb) X w 

2w 1 ^^-^^^ 

= —^b X {b ■ Vu) = - — b X [2w\\uj{ilj)b x Vz] , 

S S 

where we have used {'w)j^ = w^\b, the fact that u ■ V6 = b ■ Vw, which follows from 

B{u Vb-b- Vu) = V X {B X u) = -V X [uj{%lj)Vip] = 0, (A.17) 

I In this expression w\\ is considered to be a funciton of Rs, Ss and /i^ via (145) with aU functions of 
r evaluated at Rg in accordance with (A.G). 
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and also the fact that, for any a, 

R [{VR) (V(/)) - (V0) (VR)] ■a = axVz. 
The penultimate term in (A. 9) is gyroaveraged as follows: 



(A. 18) 



1 



[{w ■ Vb) X w] 



R 



1 



w 



wjb X {b -Vb) + ^ [V X b + {b - Vb) X b]} , (A.19) 



where we have used (A. 5). The first and second terms here are, respectively, the 
curvature drift and the so-called Bafios drift. We can use the vector identity 



b Vb = -bx (y X b) 



to write the Banos drift as 



[V X b + bb- iy X b) -V X b] 



[b ■ (V X b)] b. 



(A.20) 



(A.21) 



2 ^ 2^"^^ 
Thus, this drift is purely along the field line and one order smaller than parallel 
streaming. Therefore, as parallel derivatives of fluctuating quantities are small and 
b- V-Fqs = 0, this drift will not appear in the second-order gyrokinetic equation. Finally, 
gyroaverging the last term of (A. 9) using (A. 5) we obtain 



1 

a 



[{bxw)w-V In B] 



w1 



I R 

which is the V-B drift. 

Assembling these results, we arrive at 

,2 



21], 



b X VlnS, 



(A.22) 



b(R. 



R 



W 



^|| + ^&-(Vxb) 



(A.23) 

+ {V^)R + VBs + 0{eWj, 

where and Vbs are defined in (110) and (111), respectively (all the drifts associated 
with the mean magnetic, electric and velocity fields have been combined into Vbs). 



Appendix A. 3. Derivation of fig 

The magnetic moment Hs is defined by (66). Taking the time derivative along a particle 
trajectory, we find 

1^8 = ^ {-wlw ■V\nB + 2[v-V{bxw)]-{bxw) + 2w± ■ w} , (A.24) 

where we have used w'j_ = {b x w)-{b x w) and u-VB = (axisymmetry) . Substituting 
(A. 8) for w and rearranging terms gives 



fig = — figW ■ V In i? 



B 



-w\\ {v ■ Vb) ■ w± 
1 



(A.25) 



rris 
B 
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— yUs(if ■ V In B)j 



wn{{v ■ Vb) ■ w±)^ 



^-IIU-^ v., ^^^^^^ 

= -HsWw {b-VlnB + V -b) = -/i^wyV -3 = 0, 

where we have used (A. 4) and (A. 5). Proceeding to the third term of (A. 25), we have 
{w± ■ Vipo)ji = 0, and, substituting x = — i> ■ 6A/c, we find 



w± -Vl^ x)r + -{{w± ■ V6A^) ■ w^_) 



R 



{w± ■ V5v?)k {w± ■ V5A)^ • [wiib + u) =0, 



(A.27) 



{w±w^)j^ : Vw = 0, 



where we have used (A. 4) on each term on the right hand side. Finally, gyroaveraging 
the last term in (A. 25) we find 

(-(^.Vn)-^^)^=- 

where we have used (A. 5) and then (A. 15). 
Combining all these results, we have 

and so fig is conserved to second order in e. 



(A.28) 



(A.29) 



Appendix A. 4- Derivation of is 

The energy variable Es is defined by (65). Taking the time derivative along a particle 
trajectory, we find 

d \ (i$ 

Es = rrisV -v + Zsci — + t; ■ V [<l>(^) + <^o] - Z^e 



dt 



dip 



d 



Zscv ■{E + 6E) + Zse[— + v-V] [$(V^) + y^o] 



(A.30) 



Z.^e 



- + — i?^(..V^)^^ 



V^: + 0(e3fi,T,), 



where we have used (155) and (156) for v = as + 6as, and (69) to expand ips around ip 
in the argument of c?$/ dip. In the above expression. 



cjvj in c 

^P* = ^ + vViP+ -^R' (a, + 6as) ■ V<P ^ 
ot Z.e Z,e 



nisC ^2 



R'da^ ■ V(/>, (A.31) 



where we have used (155), axisymmetry, the axisymmetric form of the mean magnetic 
field (31) and ip = R^A ■ V0. Expanding the electric field in terms of potentials and 
using dip = u{ip)/c (see (61)), we can write (A.30) as 



ZsC ( dA dbA\ ^ 
— [n^. — ^v.—yZsev.^^^ 

—m^R? 



ZgC dip 



(A.32) 
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where we have used (32) and (62) to show that u ■ dAj dt = c 9$/ dt, so it cancelled 
with the corresponding term in (A. 30). Finally, we rearrange (A. 32) into 



Z,e f d \, Z,e dA 
+ u-V]{w- 5 A) —w 



c \dt J c dt 

Zsew ■ VSip' -^iw-Vu)-SA- '-pllR^"^ . V0) Sa^ ■ V0, 

c ZgC dip 



(A.33) 



where we have used u = u){ip)R'^'V(f) and (156). 

The lowest-order contribution to is is, therefore, 

= -Zscw^ ■ VSip' + O (e^fi.T,) , (A.34) 

and so, by using (A. 4), we conclude that 

{is) J, = O (e^nsTs) . (A.35) 

Appendix A. 5. Derivation of 

The gyrophase i} is defined implicitly by (68), viz., 

w = V — u = w\\b -\- w± (cost? 62 — sin-f^ei) . (A. 36) 

Taking the time derivative of this equation along a particle orbit gives 

as + 6as — V ■ Vu = —r^b + wuv ■ Vb + _ (cosi^Ci + smi)e2) ^. (A. 37) 

dt w± dt 

Since w± (cosi? ei + sin -(9 62) = bx w±, we take the inner product of (A. 37) with b x w± 
to find 

= — \- (as + 6as-v- Vu - w\\v ■ Vb) ■ (b x w±) . (A.38) 

Expanding according to (155) and taking only the leading-order contribution, we 
find 

^ = Vis + 0{ens). (A.39) 
Appendix A. 6. Derivation of (108) 

To derive (108) from (107), we first note that, from (A. 23), 

■ 1^ (^os + Fis + hs) = w\\b-— {Fis + hs) 

Q (A.40) 

+ n{Rs) ■ ^ + + ^Ds) ■ ^ {Fos + hs) 

because b- dF^s/ dRs = 0, and Fqs and Fis are axisymmetric. Note that the Banos drift 
(second term in (A. 23)) is negligible, as explained in Appendix A. 2. 

Now we must calculate the gyroaverages in the right-hand side of (107): 

Ts \dt[ Ts 



iss)^'-^ + ( ^ ( ^Fos ] ) . (A.41) 



R 
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This quantity is 0{e'^QsFs)- Tackling the Bohzmann response (the second term) first, 
we find 



(A.42) 

\ - \ T. -./^T^ + Ol/a.F„.), 

where we have used (105) for Fqs and the fact that OFos/ dt = O^e^Qs^os)- Using the 
lowest-order expressions (A. 10) for Rs and (A. 34) for is, and inferring from (A. 4) that 
Z^e'^{w± ■ VSip'^)j^ = 0{e^flsTg), we conclude that the second and third terms in (A.42) 
are 0{e^QsFos) and so can be neglected. Thus, we are left with 



Factoring out Fqs/Ts in (A. 41) and using (A. 33) to express ig, we find, therefore, that 
calculating (A. 41) reduces to calculating the following gyroaverage: 



is + ZsC + v - V^ 5(p' 



ZsC dip 

where x = ^V'' ~ " SA/c. 

Tackling the last term in (A. 44) first, we use the second line of (156) for 6as to find 
(neglecting O^e^Vth^Qs) contributions and so keeping only the first two terms) 

du du ^^-^'^ 

= mscR''— iycj)) ■ {vVx)r ■ V0 + (V0) ■ {w {w^ ■ V<5A))^ ■ V0. 

In the first of these terms, we split v = (^u + wi^b) +w± and use (V0)_|_ = bx ViJj/BR'^ 
(which follows from the axisymmetric form of the magnetic field, (31)); in the third term, 
we use (A. 3), integrate by parts with respect to {}, and again use (V0)_|_ = bx Vip/BR"^. 
The result is 

Tfl'^ C dill 

. y , (A.46) 

_ r|c^,d^ ^^^^ _ ^ ^^^^ _ ^ Z_^^,d^ ^^^^ ^^^^^^^ 

where (K^)j^ is defined by (110) and we have abbreviated 

duj 

{w\ib + u) ■ Vu + (Vw) ■ + u) 



W = [R^ {wiib + u) ■ V0] ^Vijj 



dip ^ (A.47) 
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We wish to write the second term on the right-hand side of (A. 46) in terms of Vw. 
Using (A. 15), we observe that 

{{b X w^) . (Vn) ■ Vx)r = R'^ m) ■ {{b X w^) V^x)r ■ V0 

+R{ib X w^) ■ [(Vi?) (V0) - (V0) (Vi?)] ■ Vxx)h- 
Rearranging this equation and using (A. 18), we find 

= {{b X w^) ■ (Vu) ■ Vx)r - {{b X w^) ■ (VxX X Vz))^ ^^'"^^^ 

= {{b X w^) ■ (Wu) ■ Vx)r + {b ■ Vz) {w^ ■ \/x)r. 

In Appendix A. 6.1, we show that {w± ■ Vx)r = and so the second term on the 
right-hand side of (A. 49) vanishes. Substituting the remainder into (A. 46), we obtain 

(A.50) 



Collecting our results, we substitute (A.50) back into (A. 44) 



is + ZsC 



d_ 
Ft 
d_ 
Ft' 



+ v 



■ V ) 6ip' 



R 



u(Rs) ■ V 



' R 



ZsC 



, dA 



— Lw-{yu)-8A) 

c 



R 



(A.51) 



-^s(y^)R ■ W + ^R'^ (VV^) ■ {w^5A)^ ■ (V0) 



where we have the fact that r — Rg = b x w/Qg to absorb the second term in the 
right-hand side of (A.50) into the first term in the right-hand side of (A.51) (this is 
valid up to corrections of order 0{e^QsTs), which we neglect). 

Turning now to the third term in (A.51), we use (A. 15) and (A. 18) to find 

{wiVu)-5A) """^"^ 



R - R'^ W) ■ {^±SA)^ ■ (V0) + u;{^){{w^ X 6A^) ■ Vz)^ 
+ {b ■ Vu) ■ (w||(5A_L - 6Ai\Wjl)j^. 



In Appendix A. 6.1, we show that the second term on the right-hand side of (A. 52) 
vanishes (see (A. 58) with a = Vz). Substituting the remainder of (A. 52) into (A.51), 
and noting the cancellation between the first term in (A.51) and the last term in (A. 52), 
we obtain 



is + Zst 



d_ 
Ft 



Z^e 



ZsC 



R 



d_ 

Ft 



+ uiRs) ■ V 



ZsC dA , 
R ^ (A.53) 



m,(Vx)^ -W + — {b- Vu) ■ {5Aiiw±_ - wn^A^) 



R' 
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It merely remains to write the first term in this expression in terms of {x)r- This 
is done in Appendix A. 6. 2. We find that 



(A.54) 



Substituting (A.54) into (A. 53), using the second expression in (A. 47) to expand W, 
and using (A. 43), we obtain our final result: 



/ ■ \ Us . 



\dt\ 


T 


TUsFos 


Iw 


T 

s 


B 



Os 



R 



d_ 
dt 



+ uiRs 



d 



OR.. 



{x)r 
dA 



(A.55) 



Appendix A.6.1. Gyroaverages involving 6A±: As V_l ■ SA_i_ = (the Coulomb gauge 
condition to lowest order in e), we can always find a scalar function C(^)^) such that 

6A± = bx V±C. (A.56) 

This expression allows us to gyroaverage explicitly various quantities involving 6A±. 
Firstly, for any vector a, 

(iu_L X 6Aj_) ■a = b a {w^_ ■ VC) • (A.57) 

Gyroaveraging both sides of this equation, we find that the right-hand side vanishes 
upon use of (A. 4). Thus, to lowest order in e, 

{{w^_ X SA^) ■ a)j^ = 0, (A.58) 

for any vector a. 

We now use (A.58) to evaluate another gyroaverage involving SA±: again using 
(A. 4), we find for x, defined in (109), that 



R 



d 



X 

The last equality follows from (A.58). 



{6Aj_ ■ {b X w^))j, = 0. 



(A.59) 



Appendix A. 6. 2. Derivation of (A.54): We start by transforming the derivative 
u{Rs) ■ Vx into a derivative with respect to r at constant Sg, Us and 'd. As ZsCU ■ Vx ~ 
0{eVLsTs), we keep track of first-order corrections. § We begin from 

2n ^ (^-60) 



+ u 



§ This is where we disagree with the derivation in [ ' ], where the difFerence between u{Rs) and u{r) is 
correctly retained, but the difFerence between u{Rg) ■ dxl dRg and u{Rs) ■ Vx is incorrectly neglected. 
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where, as the second term is e smaller than the first, we are able to drop the distinction 
between r and Rg in the argument of u. Using axisymmetry and the definitions (65), 
(66), and (68) of Eg, fis, and {}, we have (correct to lowest order in e) 

u ■ = 0, 

u ■ V/is 

u ■ 



B 



w\\ {u ■ Vb) ■ w±, 



w 



u ■ Vei) ■ 62 H ^ (■" • Vb) • (6 x w±] 



(A.61) 
(A.62) 

(A.63) 



Differentiating (109), we find 



d 



d_ 



B 



X 



X 



-5A 



1 



cragW 
1 
c 



5A^ ■ w±, 



(A.64) 
(A.65) 



Inserting these results into (A. 60) and gyroaveraging at constant Rg, we find 



{u{Rg) ■ Vx) 



R 



u(Rg) ■ V 



X 



R 



Wii 



{u ■ Vb) ■ {w±w± ■ 6A± + b X w± {b x w±) ■ 6A±) 



- {u ■ Vb) ■ {wjl6A\\). 



R 



(A.66) 



' R 



W 



where we have used (A. 17) and the fact that w±w± + {b x w±) [b x w±) = w'j_ (I — bb). 
Since Rg = r — b x w/ Qg, the spatial derivative with respect to r at constant Eg, fig 
and '(9 is converted to a derivative with respect to Rg as follows: 

„ , dx 



u{Rg) -Vl 



R 



u(Rg 



OR., 



s I R 



= u{Rs 

The last term in (A. 67) vanishes: 
u 



d{x) 



R 



dRg 



u 



V 



(b X w] 



(A.67) 



R 



[VI 



w± 

1 



(b X w) 



R 



{[{u ■ Vci) cos i}+{u- Vca) sin i}] ■ V±x)r 
{u ■ Vei) ■ e2{w^_ ■ Vx)r = 0, 



(A.68) 



where we have used the definition (68) of w, the identity (Vei) ■ 62 = — (Ve2) ■ ei, and 
finally (A. 59). Substituting the remainder of (A.67) into (A.66), we arrive at the final 
resuh (A. 54). 
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Appendix B. Flows within a Flux Surface and Moments of Fg 



Throughout this paper we need to take moments of the mean distribution function Fg. 
We need these moments in Section 7.2 to calculate the lowest-order mean current, in 
Section 8.3 to evaluate the mean flow and thence the mean Ohmic heating (198) and in 
Appendix C to prove that certain parts of the radial transport are negligible. 

In order to integrate Fg over w we need to express it in terms of r and w rather 
than Rs, Es and /i^. Assembling together the form (103) of Fg and the expressions (113), 
(114) and (118) for Fu in terms of Fi""'^ and F^*,, we flnd 



dl' B 



{E-B )^ IdA 



Fos{i){Rs),es 



+ ^ + u) ■ (V0) ^ + Fi-)(i?„ eg, fig, a 



OF 



(B.l) 



Zge ' ' ''"^^ dijj 
where the penultimate term is a convenient form for F^^ equivalent to (114). We now 
expand ip{Rg) around ip = V'(^) the flrst argument of Fog to flnd 



Fog{ij{Rg),eg) = Fog{ij,eg) + (i^, - r) 



cm,. 



(V^)^ + 0(e^F.) 
dFog 



(B.2) 



FogitP,eg) + ^R^w^-\/(l)) „^ 
ZgC oip 



Substituting this into (B.l) and combining the second term in this equation with the 
penultimate term in (B.l), we flnd 



1 - 



Zge 



dl' B 



(E-B )^ IdA 



Fog{tp,eg) 



+ "^R' {v ■ V0) ^ + Fr)(r, eg, fig, a) + 0(6^^, 



(9Fr 



(B.3) 



ZgC ^ ' dip 

where we have neglected the difference between Rg and r in the argument of f]"'^'' and 
abbreviated v = w + u. We now expand $('?/'*) around in the deflnition (65) of 

eg to obtain 
1 



eg = -mgw" - \ngU?{i,)R^ + ZgC^^ - ^i?^ [v ■ V0)' ^ + Oie^Tg 



2ng-^ ^^-^^y d^ 



where the lowest-order energy is 
1 



(B.4) 



mw + Ss = + 0{tTg 



(B.5) 



which consists of the kinetic energy of the particle motion in the rotating frame and the 
zeroth-order potential energy of the particle (centrifugal plus electrostatic) 

1 



-ITLgUJ {lp)R + ZgCifQ. 



(B.6) 
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Substituting (B.4) into (B.3) and using dFos/ dss = —Fqs/Ts gives 



Te J \ {B')^ c dt 



+ f^i?^ {V . V0) 



dip 



+ 



dtp 



where (see (105) and (96)) 



Fos{ip,eso) 



exp 



rrisW 
~2T7 



exp 



ITLsW 
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(B.7) 



(B.8) 



In the next two sections, we will use (B.7) and (B.8) to explicitly evaluate moments of 
Fs correct to first order. 



Appendix B.l. Flows within a Flux Surface 

To calculate the first-order mean flow, we multiply (B.7) by w and integrate over all 
velocities to obtain 



n,,L/.c 



dwwF, 



W {v ■ V0) 



dip \% 2) dip 



(B.9) 



nigB 



{v ■ V<pf ^Fos{iP,eso) + Fi-)(r,e,o,/is,a)| + 0(e^ 



nsVth, 



Using (B.8) for Fo,{ip,e so), we can perform the first two integrals in (B.9) explicitly: 



fd^ww^R'' («-V0) 
J Z^e 



d\nN, 



dip 



+ 



3\ c/lnTo 



dip 



FM^e 



Us 



dlnN, 



dip 



+ (S, + Ts 



dlnTs 

dip 



(B.IO) 



V0 



and 
d^ww 



rUsB 4 2 duj 

—R (v . V0) -Fo. 



-^mscR'^uj{ip)^V(p. 
ZgC dip 



(B.ll) 



Substituting these into (B.9) we find the expression for the first-order flow: 



risUs 



ZsC 

Zse 



-cR' 



Ts—r, ^ (-S + Ts) —r, \- ^sR uj{ip)-rT 

dip dip dip 



V<P+—KsB, (B.12) 
Z,e 



where we have used the fact that Jd^wwj_Fs^'^'' 
gyrophase- indep endent . 



(B.13) 

to lowest order because f]"'^-' is 
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Let us now prove that Ks is a flux function. Integrating tfie drift-kinetic equation 
(119), which determines Fs^^\ over all velocities gives 

d'wwiib-^^ = 0, (B.14) 



s 

n{nc) 



0. 



where we have used the fact that fd^ww\iFos = and fd^w C Fqs + F^^ + F< 
By exactly the same method as used in (83) to derive the lowest-order if-theorem, we 
can rewrite this constraint as 

= 2nya [ ^^B-VF]-) (B.15) 
J ml 

= ■ VK, = 0. 

Thus, Kg = Ksi^p). This immediatly implies that the first-order mean velocity is 
divergence-free (the divergence of the first term of (B.12) vanishes by axisymmetry) , 
i.e., 

V-(n,t7,) = 0(eV^^,). (B.16) 

This result will be useful in several places in Appendix C. 

Finally, multiplying (B.12) by Z^e and summing over species we find the mean 
current 

(ilnT, 



ZsCipo - ]^msUj'^{ip)R^ + Ts 



+ cR^ (^n^msR'^ u{ij)^V(f) + K{^)B 



dip 



V0 



(B.17) 



where 



K{iP) = '£k,{iP) = ld'ww\\F^-^\ (B.18) 

s s 

This is the result (126) of Section 7.2. 
Appendix B.2. Moments of Fs 

An important corollary of (B.12) is that the first-order flows stay within the flux surface: 

nJJ,-ViP = 0{e^nsV,^^ViP\), (B.19) 

which is consistent with our orderings for the timescale on which particles are 
transported. We now derive a general result that will show that toroidal momentum 
transport and heat transport are also slow across flux surfaces. We could calculate the 
lowest-order viscous stresses and heat fluxes by explicitly taking moments of (B.7) as we 
did for the flows in Appendix B.l, but we do not require these explicit expressions, only 
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the fact that their cross-flux-surface components are small. This result will be useful in 
Appendix C. 

First we prove that, for any gyrophase-in dependent single- valued function g and 
any positive integer n, 

jd^wg{r,w\\,w±,a) [R^w ■ V(/))"iu± ■ V^A = 0, (B.20) 

Indeed, using the identity (169), we find 



j d^wg{r, w\\,w±, a) [R^w ■ V(p)" w± ■ Vip 

d'^wgir^ww^w^^a) {R^w ■ V^Y — {R^w ■ V4>) 



(B.21) 



B 



n 



1 J ^'"^^ 



g{r,wi\,w±,a) {R^w ■ V0) 



n+l 



because single-valued functions must be periodic in 
We can now prove that 

Jd^wg{r,w\\,w^,a)F, {R^w ■ V0)"to^ ■ = 0(eV^7i?Xh?lWl)- (B-22) 

This follows from expanding Fg via (B.7) and applying (B.20) term by term (since w-'V4> 
differs from v ■ V0 by a gyrophase-independent function, we can replace the former by 
the latter in the statement of (B.20) where necessary). 
The radial heat flux, from (191), is 



(B.23) 



Letting g = (1/2) rrisw'^ and n = 0, we can apply (B.22) to find 

Qs-Vij = 0(eVT,t;thjVV'|). (B.24) 

We can show in a similar fashion that the radial transport of toroidal angular momentum 
(see (174)) is also small: 

(Vip) ■ Us ■ (V0) R^ = 0(eVT,i?| Wl). (B.25) 



Appendix C. Derivations for Section 8 

In this Appendix, we derive the explicit expressions (170), (184), and (196) for the mean 
fluxes in the transport equations of Section 8. We also provide detailed derivations of 
other results of that section that were stated there without proof. 

Equations (168), (182), and (C.45) express the radial fluxes as moments of dF^/ d'd. 
This is calculated in terms of known quantities via (160): 



{u + w)- VFs + 

d6f, 



Zse / ldA\ , , _ ■ 



dFs^ 
dw 



- ( 6a, 



dw 



(C.l) 



+ C[Fs] + 0{e''nsfs 



turb 
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In what follows we will use this equation to find the explicit form of the fiuxes in terms 
of hg, Fqs and Fs^^\ Note that we will use v as a shorthand for u + w, but derivatives 
are always taken at constant w. 



Appendix C.l. The Particle Flux 

The radial particle fiux is given by (168). We multiply (C.l) by (mgc/ Zgc) R^v ■ V0 
and integrate over velocities to obtain 

risUs ■ Vip = V ■ [{r\s + msUsUsU + msUsuUs + msUsUu) ■ (V</)) R^] 



(C.2) 



/ turb 



where we have used (B.16) and the identity (176) to write the first term as a full 
divergence, used ip = R^A ■ V0 to simplify the second term and used u ■ Vipo = 
(axisymmetry) . The viscous stress tensor fls is defined by (174) and 

d^wm,vC[Fs] (C.3) 

is the collisional friction force. Interpreting (C.2), we see that it relates the radial 
component of the particle flow velocity to the toroidal forces, so we can view the particle 
flux as composed of effective second-order drifts, driven by these forces. However, these 
drifts will only cause net transport of particles if they do not vanish when averaged 
over a flux surface, as indeed is obvious from the fact that we were able to average the 
particle transport equation over a flux surface. 

So we average (C.2) over a flux surface. The flrst term becomes 

^ (V ■ [(□, + m.nsUsU + msUsuUs) ■ (V0) R^] ) , = 

eld, 1 (C.4) 

Y~eV'd^V ^^^"^^ ■ ■ ^^'^^ + nisnMm'Us ■ V^>^ , 

where we have used (39) to express the flux-surface average of this divergence (recall 
also that u = u{^)R^V(f)). From (B.12), msnsUj{ip)R'^Us-Vip = 0(eVT,i?| VV'I). From 
(B.25), the same is true for the off-diagonal viscous stress. Thus, the right-hand side of 
(C.4) is 0(e^nsUthJ VV'I) and so the flrst term in the flux-surface average of (C.2) can 
be dropped. The second term does not simplify upon averaging. In Appendix D.l, we 
show that we can write the third term entirely in terms of Fqs and (see (D.14)) : 



{{F, . V0) R')^ = i^jd'w . V^) C[F, 

+ (^Jd^wF^-^^Vj,, ■ V^^^ - {n,)^m ^^^^^ 



OsJ ^ 

E-B)^ 



(C.5) 
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The fourth and final (turbulent) term can be written as follows 

\J I turb/ ^ 



73 



turb / i/) 



(C.6) 



turb/ ill 

where we have split 6fs into the Boltzmann response and the gyrokinetic distribution 
hs according to (104). Substituting from (156) for 6as, the first term in (C.6) becomes 



jd^wR^ (5as- V0) 
(V0) ■ (^Jd^w 



turb 



ZgCC 
T 

Z^ec 



V {S(p' 6 A ■ w 



{w^ ■ V) 6 A 



(C.7) 



turb 



where we have used the fact that jd^wwFos = to lowest order. The term in (C.6) 
containing hg can be manipulated as follows (again using (156) for Sa^) 



jd^^'I^R\{6a,.V<P)h,), 



turb / tp 



d^wR^ (V0) 



Vx + - (w^ ■ V) SA 



^turb^W) 



r I turb ' i/i 



(C. 



d:'w{hy^)^-^i, 



turb / i/i 



where we have used the axisymmetric form of the magnetic field (31) to write 
R? (V0) ■ Vx = (b X Vx) ■ (V'?/')/ B. The term involving bA vanished because 

[w^ . V) = ((V ■ [wMA)):)^^^,^ - (,bA{[w^ . V) 



dh 



0{en,hsSA), 



(C.9) 



'■''"II '^i/ r/ turb 

where we have used the fact that 6 A is independent of gyrophase at constant r and also 
used the identity (A. 3) in conjunction with the fact that hs is independent of gyrophase 
at constant Rg. 

Combining the above results we can write the fiux-surface average of (C.2) as 



— < n 



+ ( d-^w 



w X b 



VV- C[F, 



Os 



+ ( d'wFi^^Wns ■ ) - (n,)^/(^) 



(E-B) 



(C.IO) 



+ (^Qd'w{hsV^)^-Vij 



turb/ ip 
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The first four terms are known from classical and neoclassical theory of collisional 
transport (compare with equation 37 of [ ]) and the final term is the turbulent particle 
flux. After substitution of (C.IO) into (167), the terms containing dip/ dt cancel and we 
obtain (170). 



Appendix C.2. Amhipolarity of The Particle Flux 

It is well known that in axisymmetric tokamaks the radial particle fluxes are ambipolar, 
i.e. there is no net radial current [85, 86, 75]. This follows immediatly from the Vip 
component of the mean Ampere's Law (123) averaged over a flux surface: 

0' ■ VV^)^ = |- ((V X S) ■ ViP)^ = |- (V ■ (B X VV^))^ = 0, (C.ll) 

where we have used the axisymmetric form of the magnetic fleld, (31). 

Since j = ^s^i^sUs, to prove that our fluxes are consistent with the constraint 
(C.ll), we multiply (C.2) by ZgC and sum over species (neglecting the flrst term because 
of (C.4)): 

{j ■ vv^)^ = E (^«^) + E (^^ ■ R')^ 

^ " (C.12) 



dt I 

+ Y,^sc(^(^jd''wR^6a,-V(P)6fs'j ^ 



The flrst term in this equation vanishes due to quasineutrality of the plasma and the 
second term also vanishes because collisions conserve momentum and hence the species- 
summed friction force must vanish. The third term is 



„ \\J I turb/ X 



(C.13) 

Z^ec {R'iSnJE ■ V0),,,,>^ + {R'{Sj x 5B),^^^ ■ V0>^ , 

s 

where we have used the flrst line of the expression (156) for Sag. The flrst term in 
(C.13) vanishes by quasineutrality and so the only contribution to the radial current is 
the fluctuating Lorentz torque. But this can be written as the divergence of the Maxwell 
stress: 



R^ Q^^^l - SB6B^ ■ V0 



, (C.14) 

turb / ■(/) 

whereupon we can interchange the divergence with the average over the fluctuations 
and the right-hand side becomes ©(e^en^fthj- 

Thus, the flux-surface-averaged total radial current vanishes to second order, 
(j ■ V'?/')^ ~ 0(e^ensfthj VV'I). To the next order, the ambipolarity condition becomes 
equivalent to the conservation of toroidal angular momentum: for a more detailed 
discussion of how these two effects are connected see [751. 
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Appendix C.3. Toroidal Angular Momentum Flux 

The part of the radial angular momentum flux that we need to calculate is given by 
(182). Multiplying {C A) by {m'^c/ 2 Zsc) {R'^v ■ V0) and integrating over w gives 

(Vip) ■ n, • (V0) R^ + msn,uj{ip)R^Us ■ V^) = 



2^ 

dt 



d^w^- [R^v ■ vFs - msnsOj{i))R 
d'w^^ {R'v . V0)' C[F,] + Qd'w'^R' (v . V0) (5a, ■ V0) 5/. 



(C.15) 



turb 



where we have used axisymmetry and integrated by parts where opportune. 

Completely analogously to the calculation of the particle flux in Appendix C.l, we 
take the flux-surface average of (C.15). The first term becomes 

,2 



V- 



Jd'w^jR'vV<j^)'vF, 



V 



d'^w- 



2Z,e 



(C.16) 



where we have used (39) to simplify the flux-surface average of the divergence. From 
(B.22), we can see that we can neglect this term as it is 0{e^nsTsR\Vip\) . The second 
term in (C.15) has to do with the motion of the flux surface and will be absorbed into 
the left-hand side via the definition (167) of Tg. We split the third (coUisional) term as 
follows 

' Jd'w'^iR'vVcl))' C[F; 
Ze^ -- o (R'-^-'^^y C[Fg]} (C.17) 



c 



jd''w'^{R^wW<py C[F, 
jd^wmlR^Lo{ip) {v ■ V0) C[F, 



+ 



where we have used the fact that the full Landau collision operator conserves particles 
and so jd^w [R^u ■ V0)^ = 0. The integrals in the right hand side of (C.17) are 

calculated in Appendix D: using (D.17) and (D.12), we find 



c 
Z,e 



jd^w^ml {R^v ■ V0)^ C[Fs 
+ msOj{4j) (r^ 



w X b 



C[F, 



Os 



P^2 A^)^ii , pjBRMi^) 



(C.18) 



where ^vr,'^^''^ and ( n^'^'^^ 



are given by (185) and (186) respectively. 



Multiscale Gyrokinetics for Rotating Tokamak Plasmas 76 

Finally, we split the term involving the fluctuations as follows 

cm ^ 

* / turb 

./ Id'w''^ {v ■ V0) ((5a, ■ V0) Vi^oA (C.19) 



Vx + - {wx. ■ V) bA 

c 



(V0) K 



rl turb 

where we have used (104) for 5fs and (156) for 5as. The flrst term in (C.19), due to the 
Boltzmann response, becomes 

.2 

turb 



/C7D 

/ L3«;^^ [v . V</.) (V</)) ■ (W) S^'Fos) (C.20) 
+ //ci^iy^^ {v ■ V0) [{V5A) w-w- V5A] ■ (V0) 5ip'Fo. 



where we have used the third line of (156). The flrst term in this expression is small 
by the same logic as in (C.7) and the second term vanishes because to lowest order, 
jd^w wFqs = and Jd^wmsWwFQs = n^Tsl. As in the calculation of the particle flux, 
the term in (C.19) involving hs and Vx can be written in terms of V^: 



turb 



-( ld^wcm,R^{{vV(j)){Vx)-{'^<P)hs), 

' I d'wm^Rmw ■ V0) hsV^)^ ■ vA (C.2i; 

w / turb 

+ msLoiij)R^( [d^w{hsV^)^-Vij) 
\J 1 1 

The flnal term in (C.19) (involving bA) simplifles to 
'd^wm.R^iiv ■ V0) {yj^ ■ V(5A) ■ (V0) h,)^ 



turb 



turb 

^■''^ ^ jd'wBR'' {6 A ■ V0) (^{v ■ V0) 

[ d^wR"" (SA ■ V0) {Kw^)\ ■ V^, 

/ turb 



(C.22) 



r,tO|| ,w± 



r I turb 



where the flrst equality followed from a manipulation analogous to that in (C.9) 
and the second equality from integrating by parts with respect to "i? and using 
BR^ (dv/d^) ■ V(t> = -wi_ ■ (which follows from (31) and (68)). 
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Collecting the above results, the flux-surface-average of (C.15) becomes 
((V^) ■ n, ■ (V0) R' + ruMm'^s)^ = (vri^^% + (^r^^ 



turb / ip 



(C.24) 



- (l^jd^w^R^ {5 A ■ V0) {hgw)^ ■ Vtp 
+ (l^jd'wmsR^iw ■ V0) hsV^)^ ■ V^^ ^ 
+ msu{ij)(R''( [d^w{hsV^)^-vA ) ■ 

\ \J I turb/ V 

This will lead to the final result of Section 8.2 if we split the above expression into 
the viscous stress ti^s^'^\ as defined by (184), and the convective momentum flux 
msu{^lj) (R'Ts)^ with (R^Ts)^ given by (187): 

((v^) ■ n, ■ (V0) R' + mMm'^s)^ = {^f^^)^ + ^suj{i^) (i?'r,)^ 

\ \j ^ I turb - ^ 

Inserting this into (181) then gives (183), where the third term in (C.24) has been 
grouped together with the Maxwell stress to give the electromagnetic stresses vr^'(^\ 
defined by (188). 

Appendix C.4- Derivation of the Pressure Evolution Equation (194) 

In this Appendix, we simplify (190) and take its flux-surface average to derive (194). 

Equation (190) reads 

3d /I dA \ 

-— n^T^+V ■ Qs + ZsC i Vv?o H Qf ) ' (^^^^^ + ^s^gU ■ (Vw) ■ Us 

/ r \ ^^-25) 

+ n, : Vw - / / d^WTUsW ■ 6as6fs ) = Cf^^ + 5*^. 

/ turb 

First, denoting the potential energy of a particle S<. = Zseipo — msUj'^{'ip)R'^/2, we observe 
that 

V ■ (Q, + SsUsUs) = V-Qs + risUs ■ VEs - - Si")^ 

= V ■Qs + UsUs ■ {ZscVipo + msU ■ Vu) - R^uN)) iusUs ■ V^) ^ (C.26) 

dip 

dns 
dt 
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where we have used the continuity equation (163) to substitute for V ■ {usUs). Adding 
this to (C.25) and canceUing terms where possible gives 

i^Q-nsTs + V ■ (Q. + '^sUsUs) = - [R" (V0) ■ n, + msR'u{i^)nsUs] ■ (V^) ^ 

n.U.) - - + Ci-) + ) (C.27) 



c dt 



turb 

Before we work on simphfying the right-hand side, it is convenient to subtract 
V ■ [{usTs + Zseipo)V^], where is the velocity of the flux surfaces given by (44), 
from both sides of this equation. Then our pressure evolution equation becomes 
3 d 

-— n^T^ + V ■[Qs + '^snsUs - (n^T^ + ZsCifo) V^] = 

- [R^ (V0) • n, + msRMi^)nsUs] ■ (Vi/j) ^ 

dip 



turb 



c dt ' ' '\ dt 



V ■ [{usTs + ZsCLpo) V^] + l^j d^wnisW ■ 5as5fs 



Appendix C.4-1- Turbulent Heating. Let us first work out the turbulent contribution 
to (C.28) - the last term on the right-hand side. Using the third line of (156) to express 
5as, we can write it as follows 



d^wnisW ■ 6as6fs) = —Zgcl d^w {w ■ VSip') 6fs 

I turb 



A,e / / 3 _ , _ . . . / ,Q . „ / a 



d^wbfs {w ■ Vu) ■ 5A + d^wSfs [— + U-V \ w-6A 



c \J ' J ''\dt 

The first term of (C.29) can be rewritten by observing that 



turb 



. (C.30) 
fd^wV ■ {w6if'6fs) - fd^wSip'wwb ■ VSfs - jd'^wd^'w^ ■ V5fsj • 

J J J I turb 

Estimating the size of terms in this expression, we notice that we require 5/^ correct to 
0(e^/s) only in the last term, as the first term is an exact divergence and the second 
contains a small parallel derivative so in them we only need the first-order part of bjs 
given by (104). We evaluate the last term by using (162): multiplying it by and 
integrating over all velocities, we find 

jd^w6^'w± ■ V6fs = - jd^w6<^' (^^ + u-V^ 6fs - jd^w5ip'w\\b ■ V6fs, (C.31) 
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where on the right-hand side we now only require 6fs correct to 0{efs) to calculate the 
left-hand side to 0{e^QsnsTs/e). Substituting (C.31) into (C.30) and then (C.30) into 
(C.29) gives 



/ turb 32) 



X 



- Z.eQd'wV ■ {w5^'5fs) + Ufsw ■ (Vw) ■ 5A- 5fs{^^^+ u -V 

where we have used x = ~ ' SA/c and the fact that {d / dt + u ■ V) {x^fs)turh ~ 
0{e^flsX^fs), so it can be neglected as it is the convective derivative of a mean quantity. 
We expand Sfg in the right-hand side of (C.32) by using the decomposition (104) of Sfs 
into the Boltzmann response and the gyrokinetic distribution function: 



turb 



d^wrrisW ■ 6as6fs 

I turb 



-zJ [ d^wV ■ {6^' {whs),) + - {hsw), ■ (Vu) ■ 6 a) (C.33) 

\J ^ I turb 

where we have used Jd^wwFos = and the fact that {d/ dt + u ■ V) {6{p''^)^^j.^ = 
0(e^fis^V^'^) to eliminate the terms arising from the Boltzmann response. Next, we use 
(A. 15) and (A. 18) to work out the term in (C.33) involving Vu and find the following 
expression for the turbulent contribution to (C.28): 



/ d^wnis {w ■ 6as) 6fs ) = —ZgeV ■ ( / d^w {S(f'{whs)^ 

J I turb 

u{i:)(^jd^w{hsw), ■ {SA X V^) 



turb 



ZsC 

c 

Averaging this over a flux surface, we have 



(C.34) 



d^wnisiw ■ 6a,)6fs) ) = 

J I turb/ ?/) 

- (^(^Jd'wZ,e6v'{Kw),-Vij'^ +Pr' (C.35) 

[ d'wR^ {6 A ■ V0) {h,w), ■ ) 

J I turb/ X 



C \\J " /turb/^^V' 

where we have used (39) for the fiux-surface average of a divergence and the definition 
(204) of Pg™^. When we average (C.28) over a flux surface in Appendix C.4.6 we will 
collect the first term in (C.35) together with all other divergences to form the heat flux 
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and the final term will cancel with part of the turbulent viscous heating, as we will show 
momentarily. 



Appendix C.4-2. Viscous Heating. Next, let us handle the viscous heating - the first 
term on the right-hand side of (C.28). Flux-surface averaging this term, we have 

du 
dip 

[(Vip) ■ n, ■ (V0) R'). + ruMi^) (R'Ts)^ ^ '^'^ 



(V0) ■ n, + m,R^u{ilj)nsUs] ■ (W))^ 



+msu:{il)) I R^n 



du 



di/j\ 

' dt /^dip 



d^w{hsw)^-^SA 



du 



dip 



(C.36) 



turb 



where we have used the definition (167) of Tg, then substituted for the first two terms 
inside the brackets from (C.24), and finally introduced P™^ as defined by (197). The 
term involving 6 A is precisely the same as the last term of (C.35) with the opposite 
sign, as promised at the end of Appendix C.4.1. The last term in (C.36) will cancel a 
similar term appearing in the expression for the potential energy exchange in Appendix 
C.4.5. 



Appendix C.4-3. Ohmic Heating. We now deal with the term in (C.28) involving 
dA/dt. Using (B.12) to express nsUg and A ■ V0 = ip/R"^ (see (32)), we have 
Z.edA 



dt 



dip 
'dt 



rflniV, dliiT, 2 ( 



(C.37) 



+ K,{iP) {E + Vip)-B, 

where we have used the definition of E in terms of A and (f. Taking the fiux-surface 
average of this equation, we find that 
Z,e IdA 
~dt 



[nJJ,, 



dip 



dip 



-^+msR u{iP)- 



(C.38) 



+ p: 



Ohm 



where we have used (40) to show that {B ■ Vy?)^ = and P^^^ is defined by (198). We 
will show, momentarily, that the terms involving dip/dt cancel with identical terms in 
the expression for the compressional heating. 
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Appendix C.4-4- Compressional Heating Due to Motion of Flux Surfaces. We simplify 
the compressional heating due to V^, i.e., the V ■ [(n^Ts + Z seipojV^ij,] term in (C.28), as 
follows. Using (96) to express n^, we take the divergence explicitly to find 



' dlnN, dlnT, 2 



(C.39) 



where we have used ■ Vip = —dip/dt. We now average this over a flux surface and 
use the definition (200) of P^^^p to obtain 



(V ■ [{usTs + Zsensipo) V^]), = P; 



comp 



+ 



dijj 



dluN, 



dip 



+ (S. + Ts) 



(C.40) 



\iP) (n.V^-VP^V. 



When this is substituted back into (C.28) in Appendix C.4.6, the second line will cancel 
exactly with the similar term in (C.38) and the last line will cancel with terms from 



ppot _ shall see in the next section. 



Appendix C.4-5. Heating Due to Exchange between Potential and Thermal Energy. We 
now handle the terms on the right-hand side of (C.28) involving H^. Averaging them 
over a flux surface, we have 



„ f On 



dt 



_ 5'{") 



Zsopo - ^msU^{4j)R'^ 



Oris 

In 



Si'' 



We wish to relate this to PJ*°* defined by (201). Using (45), we can show that 

1 Oti 

2 dt 



1 



4> 

1 d 



2 V dt 



2V' dt 



V'm, (R^n,) , + —^V'-m,u?{^) I R^n 



Substituting (C.42) back into (C.41), we have 



+ 



Oris 
dt 

V'di) 



_ q(n) 



On, 



dt 



+ 



2V' dt 



V'm, (R^Tis)^ 



V'-msU^iip)(R^n 



dip\ _ 



msUj{ip) ( R^n 



du 



dip\ 

dt Ld%l) 



+ 



(C.41) 



(C.42) 



(C.43) 
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Using (201) for Pf°\ we obtain 




(C.44) 




When we substitute this into the flux-surface average of (C.28), the last term will cancel 
with the corresponding term in (C.36) and the entire second line of (C.44) will cancel 
with the third line of (C.40). The second term (C.44) will be absorbed into the definition 
(195) of the heat flux. 

Appendix C.4-6- The Flux-Surface- Averaged Pressure Evolution Equation. Taking the 
fiux-surface average of (C.28) and using (39), (45), (C.35), (C.36), (C.38), (C.40) and 
(C.44), we obtain (194). In (194), we have grouped the divergences from (C.28) together 
with the first term in (C.35) and the second term in (C.44) to get the divergence of the 
total heat flux as defined by (195). The explicit expression (196) for the heat flux is 
calculated in the next Appendix. 

Appendix C.5. Derivation of (196): The Heat Flux 

The part of the heat flux that remains to be calculated is the first term in (195), "the 
heat flux proper". By exactly the same procedure as the one employed to derive (168) 
and (182) for the particle and momentum fluxes, we can write 



In order calculate this, we multiply (C.l) by (m^c/ 2Zse) w'^R^v ■ V0, integrate over all 
velocities, and find 




(C.45) 



nisC 



[V ■ (H, + uQ,)] ■ (V0) R" 



nisC 



uj{ip)R^V ■ Qs 
T ^ 



ZsC 



ZsC 



TTLsCR^ (V0) ■ {Us + msUsUUs) ■ V(y9o - " 





R^ (V0) ■ {Us + msUsuUs) ■ {u ■ Vu) 



Vu) ■ (V0) R^ 



(C.46) 




where 




(C.47) 
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is the "kinetic-energy-weighted stress tensor" and 

= 1 j£wmlw^vC[Fs\ (C.48) 
is the "coUisional heat friction." In (C.46), we have dropped the following term: 
(fw—^R^ {v ■ V0) {w ■ Vu) ■ wFs 



ZsC dtp J 

which is 0{e^nsVth^Ts\V4'\) by virtue of (B.22) and, therefore, negligible. 

Since, from (156), it follows that (m^/ Zgc) w-6as = —w±-V6(p' + [{B / c) e^fthj 
the last term of (C.46) is 



R^l Id^w {v ■ V0) ■ 5a,5fs\ 

\J I turb 

= -m,cR^Qd^w{{v ■ V0) (w^ ■ W) hs), 

+ ^^R'jd'w (v ■ V0) . V{S^'\^^^Fo. 
= -zJ [d^w5^'{hsw^)\ ■ VV^ + 0(eVT,t;thjVV^|), 

\>'' / turb 



(C.50) 



where we have used the decomposition (104) of 6fs and the second equality is derived 
in the same way as (C.22). Using this result and adding the mean potential energy flux, 
Es = ZgCipo — msOj'^{ip)R^ /2 times (C.2), to (C.46), we get 



V ■ [m, (H, ■ V0) R^ + (n, ■ V0) i?^ + u{ij)R^E,n,Us + u{tl))R^Qs] 

^ \ Bil) r (C.51) 

T, + U,^ + — (G, + S,F,) ■ (V0) R^ 



2 J dt ZsC 

+ ^( [d^wesRHSa,-Vc^)6A , 

where Fg is the collisional friction force defined by (C.3). Taking the fiux-surface average 
of (C.51), we obtain 

(^{Qs + 'EsnsUs + Zse(^jd^w5ip'{hsw)}i ^ ■ V^^ = 
+ (( [d'we,{hsV^)^-vA ) , 

\\J I turb / i/) 

where we have used (B.25) and (B.22) to demonstrate that the fiux-surface average of 
a divergence (the first term on the right-hand side of (C.51)) is negligible (a similar 
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argument to the one that led to the neglect of (C.4)) and by exactly the same procedure 
as for the particle flux written the fluctuating heat flux in terms of hg and (see (C.7) 
and (C.8)). The collisional terms are calculated in Appendix D.l (see (D.15)): 

' - (V0)i?^X 



(C.53) 



Finally, we substitute (C.53) into (C.52) and then (C.52) into (195) to find the 
desired result (196). 



Appendix D. Collisional Transport 

In Appendix C, we derived explicit forms for the particle, momentum and heat fluxes. 
A term due to collisions on the mean distribution function Fg appeared in each of these 
fluxes and had to be expressed in terms of Fos, Fs^^\ and {E ■ B)^ (see (C.5), (CIS), 
and (C.53)). In this Appendix, we detail the relevant derivations. 

In these derivations, we will need to be able to replace {C[Fs])^ by {{C[Fs]) j^^. 
Let us prove that this is legitimate for a general function g. Taking the perpendicular 
spatial average of (A. 3), we find that the first term on the right-hand side vanishes to 
lowest order because it is a divergence, and we get 



+ Oieg), 



r^w\\ ,if _|_ 



which we can gyroaverage at constant r to find 

= 0{eg). 




(D.l) 



(D.2) 



Since any function g{r,w) can be written as g = {g)ji + dG/ d{}\ii^ for some G, we 
have, using (D.2), 





{{g)^)^ + 0(eg). 



If (7 is a mean quantity, then the perpendicular averages have no effect {{g) 
we obtain 

{{9)R)r = {9)r + 0{eg). 

Therefore, 

{C[Fs])^ = {{G[F,])^)^ 
to lowest order in e. 



(D.3) 

g) and 

(D.4) 
(D.5) 
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Appendix D.l. Collisional Fluxes 

There are three colhsional fluxes to calculate, the particle flux, the convective angular 
momentum flux and the heat flux resulting in the collisional contributions to (170), (187) 
and (196) or, within Appendix C, (C.5), (C.17), and (C.53); we handle the collisional 
viscous stress, the first term in (C.17), separately in Appendix D.2. As can be seen from 
(C.2), (C.17), and (C.46), all these fluxes can be written in the form 

Jd'wg{r,e,)iv.V^)R'C[F,], (D.6) 

where g = 1 for the particle flux, g = msR^uj{ip) for the convective angular momentum 
flux, and g = Sg for the heat flux. The calculation proceeds almost identically for each 
of these, and so we shall perform the calculation for the particle flux, highlighting the 
one point where the calculations diverge, and state the results for the other two fluxes 
at the end of this section. 

We start from the expression (C.3) for the collisional friction force, 

F, = Jd'wm,vC[F,] = Jd'wm,{w\\b + u){C[F,])^ + Jd'wm,{w^C[F.,])^, (D.7) 

where we have split v into its gyroaverage, w\ib + u, and its gyrophase-dependent part 
Wj_ (the term involving u is kept in order to maintain the symmetry between the three 
flux calculations despite the fact that it vanishes identically for the particle flux). The 
second term, the perpendicular friction, gives rise to the so-called "classical" collisional 
fluxes, whilst the first term gives rise to the "neoclassical" ones [ 13] . In the perpendicular 
friction, we can expand Fs via (103), (113), (114) and (118) to find 

[d'wm,{w^C[F,])^= [d'wms{w^{C[Fos]+ C[F*, + F^-'^]))^ 

^ (D.8) 
d^wms{w<^C\FQs])^, 



where we have used the fact that the linearised collision operator acting on a gyrophase- 
independent function is also gyrophase-independent (the neoclassical distribution 
functions only depend on the gyrophase via their slow dependence on Rg, which can 
be replaced by r in (D.8) to the order we require). Thus, we are left only with the 
perpendicular friction associated with the equilibrium distribution function (105). We 
remind the reader that Fqs is only a Maxwellian to lowest order, so C[-Fos] will be 
0{e'n,Fos)\\. 

Considering now the neoclassical friction force, we have 
d^wm,{w\\h + u){C[F,])^ 

(D.9) 

jd^wrris {w\\b + u) {{C[Fos + F^]) ^ + ( C [f]-)] >^>^ + 0(e=^ri,m,T;th.fis 



II Here we are assuming, as in the discussion about collisional energy exchange at the end of Section 8.3, 
that either all temperatures are equal or that collisional temperature equilibration is a slow process. 
Thus, there is no contribution to C[Fos] from terms like v^f'J (Tg — Tgi) F^s/Tg ^ eflgFos. 
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where we have again substituted for Fg via (103), (113), (114) and (118) and used (D.5). 
The coUision operators in (D.9) are exactly those operators that appear in (119), so we 
use (119) to substitute for them and find 

£wms{wiib + u){C[F,])^ 

^3 ( ^( J. ^^^^^^^ nrnsc jE-B)^ ^ \ (^-lO) 
= Jd wnis [w\\h + u) I ^11^ -^^-^11^-^^ {B^)J 

where we have used the fact that the integrand on the right-hand side is gyrophase- 
independent at constant r. 

Substituting (D.8) and (D.IO) into the expression for the colhsional particle fiux 
(the third term on the right hand side of (C.2)) and carrying out the velocity integration 
involving Fqs (the value of this integral will be different for the calculation of the heat 
fiux, but the procedure is identical), we find 



(D.li; 



where we have used the axisymmetric form of the magnetic field (31) to rewrite the 
classical fiux. 

For the convective fiux of angular momentum, g = msR'^uj{ip), and this is as far 
as we can proceed. Thus, the colhsional contribution to the convective fiux of angular 
momentum (the second term in (C.17)) is 

c 



d^wmlR^u{ip) {v ■ V0) C[Fs] 
+ mMi^)(^R^ jd^'w (^^f^-VV^^ C[F, 



(D.12) 



■ 

This expression is used in (CIS) and then in (187). 

For the particle and heat fiuxes {g = 1 and g = Eg, respectively) we can make one 
further simplification. Using Jd^w = J desd^sd^{B /m1w\\) and integrating by parts 
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jd^wln^ (V0) ■ {w\\h + u) w\\h 



OR, 



(D.13) 



wfiere we liave used tlie identity (116) obtain tlie last expression. Substituting this into 
(D.ll), we obtain the final form of the collisional particle flux 

(D.14) 



c 



(E-B) 



which is (C.5) and contains the collisional terms appearing in (170). 

As all the derivatives in (D.13) are taken at constant Ss-, the derivation presented 
for the particle flux {g = 1) can be carried out identically for the heat flux {g = Eg) (the 
fact that b ■ V [msR'^uj{il!)] 7^ is what prevented us from using (D.13) to simplify the 
convective angular momentum flux). Thus, the heat flux is 



Z.e 



{{Gs + S,F,) ■ (V0) R')^ = Y^j ^'^^^ • R^C[F, 



+ ( d^WBsF^^'^Vr^s-V^ 



(D.15) 



which is (C.53) as required. 



Appendix D.2. Collisional Viscous Stress 

We now calculate the collisional contribution to the viscous stress, the first term in 
(C.17). Splitting {R'^w ■ V(j)f into its gyroaverage and a gyrophase-dependent part, we 
find 

Z.e 



jd'w'^{R^wVct>f C[F, 



d^w'^{{{R'w^.V<py 



[B'R'-I\^)]\C[F, 



r I i/j 



+ 



C 

Zse 

c 
Zse 



d^wml^-^^R\{w^.V<P) C[F,s])^ 



{{C[Fs])^)^ 



(D.16) 
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where we have used {{R^w ■ V(f)f)^ = [P{ij)/B^] U7j + [i?^ _ P{iIj)/B^] wl/2 and (D.5). 
The first two terms in (D.16) are the classical contribution to the collisional viscous 
stress. They can be rewritten using the axisymmetric form of the magnetic field (31) 
and thus shown to be equal to /vri'^''') as given by (185). The last of the terms in (D.16) 



is the neoclassical viscous stress. In the same way as for the neoclassical particle flux, 
we use (119) to substitute for {C[Fs])j^ and integrate by parts under the fiux-surface 
average to find that it is equal to (ttI"'^^ ) as given by (186). 



Thus, the collisional viscous stress is 



c 
Z,,e 



jd'w'f {R^w . V0)^ C[F,]j^ = + , (D.17) 

where the classical and neoclassical parts are given by (185) and (186), respectively. 

Appendix E. Derivations for Sections 9 and 10 

In this Appendix we provide detailed derivations of equations that are stated without 
proof in Sections 9 and 10. 

Appendix E.l. Derivation of (210): Evolution of Thermal Energy 

Summing (194) over all species we find 

3^ d_ 
2V' di 



''P s s 



pji- ^ pOhm ^ pcomp ^ Y °* + Pf 1^ + Y {Sf^) 



(E.l) 



where we have used the fact that collisions conserve energy and so Y.S Cs^^ = 0- We 
now proceed to simplify the source terms on the right-hand side of this equation. 
Starting with the viscous- heating term, defined in (197), we find 



where we have used (183) to substitute for the momentum flux in terms of the total flux 
and its electromagnetic part. 

We handle the Ohmic-heating and compressional-heating terms together. Adding 
(C.38) to (C.40) and summing over species, we obtain 



(E.3) 



s 
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where we have used (44) for Vip and used quasineutrahty to ehminate all terms 
involving cpQ. Rearranging this equation, we have 

(E.4) 



s s ^ ' 



32 

where we have used (39) for the fiux-surface average of a divergence and used 
— {l/c)dA/dt = E + V<f. The first term on the right-hand side of (E.4) can be written 

as 

where we have used V-j = 0, the identity (39) again, ip = ^{ip) + ipQ, and (j ■ Vip)^ = 
(see (C.ll)). 

Finally, turning to the potential-energy term, we collect the second term in the 
last line of (E.4) together with the potential-energy-exchange term, and substitute the 
definition (201) to find: 



{{E + Vv) . j)^ = {E . j)^ + (V ■ i^j))^ = {E . j)^ + -—V {^oj ■ VV^)^ , (E.5) 



2V' dt 



(E.6) 



where the moment of inertia J is defined in (180) and the terms involving (po vanished 
by quasineutrahty: ZgCUs = 0, ZgeSi"^^ = 0. 

Substituting (E.2), (E.4), (E.5), and (E.6) back into (E.l), we find the desired result 
(210), where we have grouped the second term in (E.4) and the last term in (E.5) with 
the divergence of the heat flux on the left-hand side of (E.l) and also grouped the particle 
source term from (E.6) with the heat sources in (E.l). We postpone the discussion of 
the turbulent-heating term P*"^*^ until Appendix E.2, as it is particularly involved. 

Appendix E.2. Derivation of (211): Turbulent Heating 

Here we calculate the contribution to the evolution of the thermal energy due to 
fluctuations, i.e., the source term '^gPg^'^^ on the right-hand side of (210). Starting 
from the definition (204) of P*"''^ we have 

7 I I {■ \ \ ^ ' 

s ^ \\J I turb/ V 

When we come to derive the conservation law for free energy in Appendix E.4, we will 
need an expression for the right-hand side of (E.7) except without the time average that 
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is implicit in the average over fluctuations (see the definition (18) of the fiuctuation 
average). We will therefore first carry out our calculation of the right-hand side of (E.7) 
without this time average and then apply the time average at the last step to obtain 
J2s Pg^^^- As Section 9.2, we refer to the composition of the perpendicular spatial 
average with the fiux-surface average as the annulus average. 

Thus, averaged only over an annulus but not over time, the first term on the right- 
hand side of (E.7) becomes (using the perpendicular spatial average as defined by (15)) 



c\\ \dt I /./, ^\dt\ 2T. 

(6E + V V) -Sj + -(ux SB) -Sj + - (6j ■ Vu) ■ 6 A 

c c 



(E.8) 



where we have used hs = Sfs + Zgcdip' Fqs/Ts, x = ^v' '~ " ^A/c, the quasineutrality 
constraint (146), and the fact that (w-V(---)_|_)^ = 0, which follows from V ■ u = 0, u ■ 
Vijj = and the formula (39) for the fiux-surface average of a divergence. The last line in 
(E.8) is obtained by using 86 A/ dt = -c {6E + \/6ip) = -c {6E + VV) - V ■ (w ■ 5A). 

We now find an expression for {{SE + VSip') ■ We start from the 

perpendicular spatial average of Poynting's theorem for the fiuctuating fields (219): 

^ ^^^^^ +^\/-{6Ex6B)^ = -{6E-6j)^. (E.9) 



dt \ 8tt I ^ Att 

Taking the cross product of (98) with 6B, we have the following expression for the 
Poynting fiux: 

—SE x6B= —u ■ (55^1 - 6B6B) - — (VV) x ^B. (E.IO) 
An An An 

Using (E.IO), we obtain 

I + ^ ■ (^)^ - ■ ■ ^) = - + ^'^^') ■ ■ (E-ii) 

Using this to substitute for {{5E + \/6(p') ■ 6j)j_ in (E.8), we find 



X 



^6B6B :Vu+- {6j ■ Vu) ■ 6a\ \ (E.12) 

d I ZyS^'^ris \ \ I d I bB" 



^\dt\ 2T, /,/.. \dt 



<>n 
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where we have used Ampere's law to write [1/ c) {{u x SB) ■ Sj) j_ 
- (l/47r) {6B ■ {V6B) ■u)^^. Using (A. 15) and (A. 18) to express Vw, we obtain 



d^wih, f ^ + V ) X 



( VV^) ■ ( ^5B5B + -5j5A\ ■ (V0) i?^ 



An 



dip 



(E.13) 



dt 



2T. 



d /5B^ 



dt \ Stt 



Using 6j = '^gZge jd^wihgw)^ in the second term on the right-hand side of this 
equation, we arrive at the expression for the non-time-averaged turbulent heating: 

d 



c 



(VV^) ■ ( j^5B6B + -JjdA"^^ . (V0) R'^ ^ 



(E.14) 



Ft 



2T, 



d /5B^ 



dt \ Svr 



where we have moved the term involving 6A x Vz to the left-hand side to make the 
comparison to (E.7) more transparent. 

Finally, taking the time average of (E.14), we obtain from (E.7) 



turb 



(E.15) 



where we have used the definition (188) of and the fact that time derivatives of 
averaged quantities are smaller than the other terms and can be neglected. This is 
precisely (211), as required. 



Appendix E.3. Derivation of (215): A Simple Form of the Kinetic Energy Flux V'^^^ 

Starting from the first expression for V'^^'^ of (215), we wish to derive the second. Using 
(195) to express {qs)^ in the first line of (215), we obtain 



(E.16) 
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where we have used j = ZgengUs and so cancelled the last term in the first line of 
(215) with corresponding term in (195). Next, we expand vr^ot''^'' in (E.16) by using its 
definition (181) and obtain 



turb. 



(E.17) 



where we have used u = uj{ip)R^V(t) and the symmetry of n^. We now use (167) to 
express Tg in terms of ngUs and dip/ dt and find that 



turb. 



turb 



(E.18) 



where we have collected all terms proportional to ngUg and dip/ dt together and used 
the definition (214) of the total kinetic energy U . 

To simplify the terms involving the fluctuations, we need to relate u • {5B5B)^^^^ 
to the Poynting flux. Averaging (E.IO) over the fluctuations, we obtain 

""-{dE X 6B)^^^^ ■ VlP - - — ^ - W-.-/t,,b - V - — \v V .y. ; ^ w^/turb 



An 



~u ■ {6B6B),^^^ . ViP - ^((VV) X " V^- (E.19) 



(E.20) 



Using the fluctuating part of Ampere's law (148), we rewrite the last term as 
-^((W) X 5B),^^^ ■ ViP = Y,Zsel^jd^w5^'{Kw)^^ ^ ■ ViP 

- ■ (^V'^B X V^)turb- 

The second term in the right-hand side of (E.20) is order e smaller than the flrst and 
so we can neglect it. Thus, substituting (E.20) into (E.19), we flnd that the sum of the 
two fluctuation terms in (E.18) is equal to (c/ An) {6E x 6B)^^^^^ ■ Vip. This gives us 
the second line of (215), as required. 



Appendix E.4- Derivation of (222): Free Energy Balance 

In this Appendix, we derive the free-energy evolution equation (222) from the gyrokinetic 
equation (144). The free energy is deflned in (221). We start by using the expression 
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(104) for 6fs in terms of hg to rewrite W as 



w- 



+ 



53' 



2T, 



Stt 



(E.2i; 



where we have used the fluctuating quasineutrahty constraint (146). 

We now find an evolution equation for the first of the terms in (E.21). Multiplying 
the gyrokinetic equation (144) by Tghs/ Fqs, we obtain 

d_ d 
di 



u(R, 



2F, 



Os 



dRs 2Fos 



Z.eh. 



d_ 
Ft 



+ u(Rs 



d 



dR: 
B 



{x)r 



(E.22) 



du 
dip 



We now average this equation over an annulus and integrate over all velocities. 

We start by processing the left-hand side of (E.22). Using the identity (45), the 
time derivative becomes 



d^w 



dt 2F, 



Os 



d_ 
dt 



d^w 



2F 



Os 



(E.23) 



where have dropped the term arising from the motion of the flux surfaces because it is 
smaller than the rapid variation of h"^ . We have neglected the time variation of V on 
the same grounds. In Appendix E.4.2, we show that 



j d^w(^ [u{Rs) +wi\b + Vds + {V^)r\ 



OR, 2F( 



Os 



Oie^n^W). (E.24) 



Thus, the various advection terms in (E.22) do not contribute to the evolution of free 
energy. Note in particular the vanishing of the term involving (T^)^ - a reflection of 
free energy's fundamental status as a cascaded invariant with respect to the gyrokinetic 
nonlinearity. 

Consider now the right-hand side of (E.22). The first term can be interpreted as the 
source of (equivalently the perturbed entropy) due to turbulent heating. In Appendix 
E.4.3, we show that the average of this term is 



d wZ,e( h 



d_ 
dt 



+ u(Rs 



d 



d w( h 



d_ 
dt 



dR, 



R 



ZsC 

c 

Z^e 



Jd^w{hsw)^ ■ (5 A X V^) 



(E.25) 



d^w( h 



b X w 



Multiscale Gyrokinetics for Rotating Tokamak Plasmas 



94 



Finally, we turn to the second term on the right-hand side of (E.22) and calculate the 
source of due to background gradients (equivalently, free-energy injection due to 
instabilities). In Appendix E.4.4, we show that 



d'w{Ts—^ + ms 



B 



jd^we,{h,V^)^-Vij 



dip 

d\nNs 3d\nTs 
dip 2 dip 

dlnTs 

dip 

du 

dip 



(E.26) 



d^'wi h, i • ) ■ Vx 



Collecting (E.23), (E.25), and (E.26) together, we find that the annulus average of 
(E.22) integrated over all velocities is 



d_ 
dt 



d-^w 



2Ff 



Os 



[d'w^{hsC[h])^ 

J ^Os 



d 



Zse 
c 



d^wi h, { — + U-V ]x 



u{iP) (^Qd'w{h,w)^ ■ {6A X V2 



d'w{KV^)^-Vip) ) T, 



d\nN, 3dlnT, 



dip 2 dip 



(E.27) 



d'weAh.V,)^.Vip) ) ^ 



d'wR' {V(p) ■ (vhsV^)^ ■ Wip 



du 

dip ' 



where the last term in (E.25) and the last term in (E.26) have cancelled. Summing 
(E.27) over all species, using (E.14) for the first two terms on the right-hand side, and 
collecting together the time derivatives, we obtain (222). 



Appendix E.4-1- Derivation of ( 204 )■ Averaging (E.27) over the intermediate timescale 
T and using the fact that (fi')^^,]-, = {{g)±)rp, we find 

d 



jdiss 



r I turb ' i/) 



Z^e 



Uj{lP) 



d^w{hsw)^ ■ {6A X V^) 



turb/ ip 



(E.28) 



+ p: 



drive 
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where P^^^^ and P^^"'^ are defined by (205) and (206), respectively. Rearranging terms 
and using the first line of (204) (the definition of P*"^*^) gives the second hne of (204) 
immediatly. 



Appendix E.4-2. Derivation of (E.24)- Here we prove that the advection terms in 
(E.22) do not contribute to the evolution of the free energy. Starting with the mean 
rotation, we have 



Os 



s"'s 

2-Fns 



d 



dR.. 



Tshl 



2Fa 



(E.29) 



where we have used (A. 60) and (A. 61) to change the spatial derivative from one at 
constant e^, /i^ and 'd to one at constant w. Substituting for u ■ V/x^ and u ■ V'd from 
(A. 62) and (A. 63), respectively, carrying out the gyroaverage, and integrating by parts 
with respect to d where necessary,^ we find that the d/ dfig and 8/ dd terms in (E.29) 
vanish and we are left with 



d w{^{Rs 



OR, 2Ff 



Os 



d wu(R. 



d 



OR. 



2Ff 



Os 



(E.30) 



We now need to rewrite the derivative with respect to -R^ as a derivative with respect 
to r in order to interchange the derivative and the perpendicular spatial average. Using 



d 



u 



OR, 



u V + u 



dr d 



OR, OR, 



u 



(Vb) X w 



(E.31) 



and the definition (68) of w to do this, we obtain 
d^wi u(R, 



OR, 2Fr 



Os 



jd'w(u{R,)-V^'^ ^ -{u.Ve,).e,jd'w(^-V 



TMl 



2Fr 



Os 



(E.32) 



Note that, as we have been able to interchange the derivative and the average, the first 
term in (E.32) is now one power of e larger than second term, and so we can drop the 
latter. Similarly, we can drop the distinction between Rs and r in the argument of u 
and fiux-surface average (E.32) to obtain 



d^w(u{R,)-^'^ 



Os 



2Fn 



(E.33) 



where we have used V ■ it = 0, (39) for the fiux-surface average of a divergence, and 
^ Note that we are allowed to interchange the ()^ and averages. 
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Turning to rest of the advection terms in (E.22), we find that they can be rewritten 
in divergence form: 



I I r I Th'^\ \ \ (EM) 

= {V ■ [jd'w{^ h|b + Vd. + (V,)^) jj^ + Oie'QsW), 

where we have used the fact that derivatives with respect to Rs and with respect to r 
are equivalent to lowest order, used an analogous calculaiton to (83) to interchange the 
divergence and the velocity integral of wyb- V, and used V- = 0{e^fls)- Using (39), we 
see that the term involving w\\b vanishes. As the spatial derivatives of a perpendicularly 
averaged quantity are small, the remaining terms in the right-hand side of (E.34) are 
0{e'^QsW) and are thus negligible. Combining this with (E.33), we obtain (E.24) as 
required. 

Appendix E.4-3. Derivation of (E.25). Using (A. 54) to substitute for u{Rs) ■ 
9Rs), we have 



Z,e( ( Id^wi h 



u{R, 



dt ' dR., 

d_ 
dt 



c 

\ B 



c 



^(^jd''w{hs{bxw)-{Vu)-Vx)r 

b ■ Vu) ■ l^j d^w{hs {w±5Aii - wi\SAj_) )^ 



(E.35) 



where we have expanded u{Rs) about u{r) and used (D.3) to change from gyroaverages 
at constant Rs to gyroaverages at constant r. The last term in (E.35) can be written 

as 

^(^(b-Vu)- /^jd'w{K K^ij - w\\5A^))^ ^ 

= -^uj{^) (^jd^w {{hsw)^ ■ {6A X V^))^^ (E.36) 

+ ^uii:) ^(6 ■ Vz) jd'w {{hsb ■ {w^ X 6A^))^) 

where we have used (A. 15) and (A. 18) to expand b ■ Vu. Using (A. 56) to express SA±, 
we see that 

{{Kb . {W^ X 6A^))^)^ = {{KW^ . VOr)± = V ■ {{w^KOr)^ 

= 0{evth,hs6Aj_), 

where we have used (A. 4) and the fact that hs is independent of at constant Rg. 
Thus, the second term on the right-hand side of (E.36) can be neglected. Substituting 
the remaining term back into (E.35) gives (E.25) as required. 



-L / ' 
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'MV^)R)^)={(hsV^)r)^- 
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(E.38) 



Therefore, 



dip 

Jd^wT,{KV^)^-Vtlj 



din AT, 3dlnT, 



dip 2 dip 



(E.39) 



dlnTg 

where we have used (105) for F^g. Turning to the term in (E.26) involving du/ dip, we 
subtract (A. 50) from (A. 45) and rearrange terms to find 



m.ci?^^ (V0) ■ {vVx)r ■ V0 



where W is defined by (A. 47). Therefore, 



b X w 



(E.40) 



R 



— m„ 



d^w 



I{ip)wii 



B 



+ co{iP)R' 



du 
dip 



d^w{Km,{V^)^-W) 
Jd^wR^ {V(P)-{vhsV^)^-ViP 



du 
dip 



(E.41) 



+ 



nisC 



d^w{h,{b X w) ■ (Vw) ■ Vx), 



where we have used (A.47), b x Vip = BR^Vcp and (D.3). Adding (E.41) to (E.39), we 
obtain (E.26) as required. 



Appendix E.5. Derivation of (231): The Entropy Flux 

Here we derive an expression for the radial entropy fiux in terms of known quantities. 
Letting fs = Fs + 6fs in the definition (228) of this fiux and then expanding Fg = 
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Fqs + Fis + ■ ■ ■ inside the logarithm, we obtain, keeping terms up to 0(e^ 



;r'"'>, = -E(/*'«' 



FAn 



STc^h^F 



mi 
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w± ■ V^p 
d^wF, <( 1 + In 



+ Fu + F2s + 



F2 



2F 



Os 



""mi 

2txK 



2 \ 3/2- 



T, 



(E.42) 



turb 



^(^jd^wFosWx-Vij 



In deriving the last expression, we are helped by the fact that only the gyrophase- 
dependent parts of the distribution function give non-vanishing contributions - so we 
have used the gyrophase-independence of Fu (see Section 7.1), (104) and (D.3) to 
simplify the term involving 5/^, ehminated F2s by writing F^ + F2s = Fg — Fqs, and 
used (B.7) and (B.22) to ascertain that we do not need to keep O(eFs) corrections to 
Fos when calculating InFos. Expanding Fos{'ip{Rs),£s) about a local Maxwellian, as in 
(B.2), and retaining terms up to O(e^Fos), we can apply (B.20) multiple times to show 
that 



d^wFos'w± ■ Vip 



0(e\th.ri.|VV>|) 



(E.43) 



and so the final term in (E.42) can be dropped. Using (96) inside the logarithm in the 
lowest-order expression for the mean entropy density (225), we find 



2 \ 3/2- 



+ 



2T. 



ZsClpq _ 3 
T. 2 



(E.44) 



which we substitute into the last term of (E.42). Now, using (167) to write 
jd^wFgWj^ ■ Vip in terms of F^, (B.4) to expand e^, and working to the lowest non- 
zero order in e, we have 



;r(H)) 



s 
s 
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2 \ 3/2- 
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Using the definitions (191) of Qs and (165) of ngUs in the definition (195) of Qs, we see 
that we can rewrite the second and third hnes of (E.45) in terms of to arrive at 

(r"">v, = -E + 

Using (232) for the chemical potential results in (231) as required. 



2\ 3/2 



(E.46) 



Appendix E.6. Derivation of (235): Collisional Entropy Production 

In this Appendix, we derive (235), expressing (233) in terms of known quantities. In 
particular, we wish to write as many terms as possible in terms of the transport fluxes 
derived in Appendix C. To this end, we start by splitting the fluxes into collisional and 
turbulent contributions. 

Using (170) and (C.5), we can write the particle flux in terms of its collisional and 
turbulent contributions as follows 

{Ts)^ = ^({Fs-V<P)R') +(( fd'w{hsV^)^-vA ) . (E.47) 

■^sd \ \J I turb/ 

Similarly, from (196) and (C.53), we have for the heat flux 

(g.)^ = ^((a + S,F,)-(V0)i?2> +///rf3^£.(/i.V^),-V^\ \ . (E.48) 

\\J I turb/ V 

For the momentum flux, we first use (CIS) in (C.23) to obtain 
((W) ■ n, • (V0) R^ + m,uj{^)R^Y,)^ = 

J d?w^R^ {5 A ■ V0) {Kw)^ ■ 
' [d'wm,R^{{v ■ V0) h,V^)^ ■ 

W / turb/ ip 

where we have used v = u + w to group the last two lines of (C.23) together. Now, 
after we use (C.24) to rewrite the left-hand side of (E.49) in terms of iri^'^^ and F^, the 
term involving 6 A cancels and we find 



c 



turb / ip 



(E.49) 



R'^s), =^ Qd'wlml {R'v . V0)^ C[F- 



+ 



I turb/ ih 



(E.50) 



turb / i/i 

We will now proceed to show that the first term in (233) can be expressed via the 
collisional terms in (E.47), (E.48), and (E.50), and that the second term in (233) can 
be expressed via the turbulent terms in (E.47), (E.48), and (E.50), with precisely the 
same coefficients. 
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Using (B.7) to expand Fg inside the logarithm in the first term in (233), we obtain 



+ 



cnis 



mi 



d^w 



c[f; 



dlnNs 

dip 



+ 



d^w ( In 



+ 



TUsB duj 




8iT^h^Fos{ilJ,eso) 



C[F- 




(v-Vcp) C[F, 



d w 



(E.51) 



21],T, dip 

where we have used the fact that colhsions conserve particles (this eliminates the dA / dt 
term arising in (B.7)) and the fact that CfF,] ~ e^Qs^s because u ~ eQg and F^ 
is Maxwellian to lowest order. Multiplying the neoclassical kinetic equation (119) by 
Fs^^'^ / Fqs, integrating over all velocities, and flux-surface averaging, we obtain 



^d^w 




pOhr 

s 



(E.52) 



where Kg is defined in (199) and we have used a manipulation analogous to the one in 
(83) to eliminate the first term in the left-hand side of (119). Substituting (E.52) into 
(E.51), and making use of (B.8) for Fos{ip,£so) and the definitions (192) of the collisional 
energy exchange Cs^\ (C.3) of the friction force i^^, and (C.48) of the "collisional heat 
friction" Gg, we find 

^Ti^h^F, 



d^whi 



C[FP 
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dip 
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^ s 

3 (i In Tg 
" 2 dip 
c/lnT, 



Ohm 



(E.53) 



dip 



where Eg = ZgCipo — {1/ 2) nisu'^R'^ is the potential energy of a particle as in (B.6). 
We see that the moments of the collision operator in the last three lines of (E.53) are 
precisely those in the expressions for the fluxes derived above. 

We now complete the derivation by obtaining the fluctuating part of the fluxes from 
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the second term in (233). We have, by definition (205), that 



\\J ^Os /turb/ ,6 -'-s -l-s 




dlnNs 3d\nTs 
dip 2 dip 



) 



(E.54) 



+ ^(( [ d'wmsR' ( V0) ■ {vKV^)^ -Vip) ) 




dip ' 



duj 



I turb / i/j 



where we have used the second hne of (204) and then (206) to write p^^"™ in terms of 
the gradients of n^, T^, and uj. We see that, again, the turbulent fiuxes in (E.54) match 
up precisely with the turbulent fluxes in (E.47), (E.48), and (E.50). 

Substituting (E.53) and (E.54) into (233) and using (E.47), (E.48) and (E.50), we 
arrive at (235) as required. 
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